{'fraction_non_alphanumeric': 0.056914794104808035, 'fraction_numerical': 0.018250982185544724, 'mean_word_length': 4.224019980453904, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On port-Hamiltonian approximation of a nonlinear flow problem on networks

This paper deals with the systematic development of structure-preserving approximations for a class of nonlinear partial differential equations on networks. The class includes, for example, gas pipe network systems described by barotropic Euler equations. Our approach is guided throughout by energy-based modeling concepts (port-Hamiltonian formalism, theory of Legendre transformation), which provide a convenient and general line of reasoning. Under mild assumptions on the approximation, local conservation of mass, an energy bound, and the inheritance of the port-Hamiltonian structure can be shown. Our approach is not limited to conventional space discretization but also covers complexity reduction of the nonlinearities by inexact integration. Thus, it can serve as a basis for structure-preserving model reduction. Combined with an energy stable time integration, we numerically demonstrate the applicability and good stability properties of the approach using the Euler equations as an example.

1. Introduction. Structure-preserving discretization is an active research area in the last decades. By preserving or mimicking relevant geometric structures such as, e.g., conservation laws, dissipative relations, or symplecticities, unphysical solution behavior and numerical instabilities can be avoided in many cases, cf. [42,18,14,2]. The model problem we consider in this paper describes nonlinear flows on networks. It covers a hierarchy of models used to describe gas network systems, including particularly the barotropic Euler equations [52,20,21,25,19], but also p-systems [61,51] and more general symmetrizable hyperbolic systems [36,53]. Further applications are in the contexts [40,48], e.g., in the modeling of electric transmission lines. The network is assumed to be described by a directed graph. Each edge ω of the graph can be identified with an interval. Given a strictly convex smooth function h : R 2 → R and a non-negative functionr : R 2 → R, the edgewise statesz ω = [z ω 1 ; z ω 2 ] : [0, T ] × ω → R 2 are governed by ∂ t z ω 1 (t, x) = −∂ x ∇ 2 h(z ω (t, x)), ∂ t z ω 2 (t, x) = −∂ x ∇ 1 h(z ω (t, x)) −r(z ω (t, x))∇ 2 h(z ω (t, x)), x) for i = 1, 2. The expressions M(z) = ω∈E ω z ω 1 dx andH(z) = ω∈E ω h(z ω )dx, with E the set of all edges, represent the total mass and the Hamiltonian of the system. Fundamental properties of the hyperbolic model problem are that, under appropriate coupling conditions on the edgewise equations, conservation of mass and dissipation of the Hamiltonian (energy dissipation) hold up to the exchange with the boundary. Moreover, the convective terms can be related to a certain skew symmetric geometric structure. Our aim is to derive a structure-preserving approximation approach, which is applicable for both conventional space discretization and projection-based (Galerkin) model order reduction with additional complexity reduction. The latter denotes a sparse approximation of nonlinear terms (also known as hyper-reduction). Note that we use the expression model reduction as an umbrella term for model order-and complexity-reduction.
There exists a rich literature on structure-preserving discretization. The mimetic finite difference and finite element methods [6,41] are designed to mimic conservation laws and dissipative relations on the discrete level. These approaches typically involve the separate approximation of operators and geometric objects and do therefore not completely fit in the Galerkin framework, or are tailored towards specific discrete structures or problems, see e.g., [21,58,26,25] in the gas network context. Many entropy-stable methods are based on a flux-centered point of view and the so-called entropy-flux pair. This is the case for finite volume methods [42,39] or methods using summation-by-parts approximation [29,63]. These methods are particularly wellsuited in the presence of shocks. The summation-by-parts methods have the additional advantage that inexact integration can be taken into account, but they are not purely projection-based and heavily rely on the flux-interpretation, which complicates their adaption for model reduction. A contribution in this direction we are aware of is [16]. For systems in Lagrangian form, structure-preserving Galerkin approximations are studied in [13,27]. The Hamiltonian formulation is dual to the Lagrangian one [18,30], and the so-called Legendre transformation links the Hamiltonian and the Lagrangian functions of a system [59,68].
For the approximation of our model problem, we consider the port-Hamiltonian framework, a generalization of the Hamiltonian formalism, which is particularly wellsuited for network problems. It has its origins in the analysis of finite-dimensional connected systems [47,48,65,64] but has been extended to the infinite-dimensional setting of partial differential equations [7,31] and has been systematically generalized to constrained dynamical systems [9,64]. The lumped port-Hamiltonian approximation of infinite-dimensional systems has also gained interest recently. Most contributions in this direction either focus on linear systems [34,67] or are tailored to very specific discrete structures such as mimetic discretizations [29,28,56]. There exist a few works on structure-preserving model reduction of nonlinear port-Hamiltonian systems, most of them consider the finite-dimensional case, see [17,2,3,57]. A comprehensive literature is available for standard model reduction. cf., e.g., [8,10,27,13].
In this paper, we propose a conforming Galerkin ansatz, which respects the Hamiltonian and geometric structure of our nonlinear flow problem. On the one hand, we make use of compatibility conditions on the ansatz spaces, which also play a fundamental role in mimetic finite element methods and symplectic model reduction [57,2]. On the other hand, a certain variable transformation related to the Hamiltonian of the system and its systematic analysis using the theory on the partial Legendre transformation is crucial in our approach. To the best of the authors' knowledge, the latter has not yet been used in this extend in the context. The use of structured variable transformations widens the range of formulations, for which structure-preserving Galerkin approximations can be derived with purely variational arguments. Moreover, complexity reduction of nonlinear terms by inexact quadrature can be included without much difficulty, similar to other Galerkin-based approaches [22,27]. We present a complete modeling work flow consisting of the following steps: Port-Hamiltonian modeling of our flow problem on networks; analytical investigation of a variable transformation induced by the partial Legendre transformation; approximation in space by Galerkin projection and quadrature-based complexity reduction (hyper-reduction); energy-stable time discretization.
The structure of this paper is as follows: The underlying energy-based modeling concepts, i.e., the port-Hamiltonian formalism and the partial Legendre transformation, are shortly introduced in Section 2. In Section 3 we present our model problem together with an appropriate parametrization of the solution and a variational principle. The latter provides the basis for the structure-preserving approximations derived in Section 4. The approximations inherit, among others, port-Hamiltonian structure, which is revealed by their structured coordinate representations in Section 5. An energy-stable time discretization is presented in Section 6 and the applicability of our approach is numerically demonstrated for the barotropic Euler equations in two test cases in Section 7. Our conclusion and outlook give a summary and point towards possible directions for future research. The appendix extends our approach to systems with an additional dissipation term as well as with edge weighting.
Notation. In this paper matrices, vectors and scalars are indicated by capital boldfaced, small boldfaced and normal letters, respectively, whereby vector-expressions always refer to column vectors. Given two vectors a ∈ R n and b ∈ R m , we write [a; b] ∈ R n+m for the concatenation to a new (column) vector. For m = n the Euclidean scalar product is denoted by a · b. We distinguish between functionvalued vector spaces (e.g., L 2 , V) and real-valued spaces and sets (e.g., R N , S) in the typesetting.
Given a scalar field h : R n → R and a partitioning of the argument into vectorcomponents z = [z 1 ; z 2 ] with z i ∈ R ni , n 1 + n 2 = n, we often structure the gradient ∇h and Hessian ∇ 2 h in sub-blocks, i.e., The solutions of the partial differential equations in this paper depend on time t and space x, e.g., z(t, x) ∈ R. When it is convenient, we interpret them as functions in time with values in a function space, e.g., z(t) ∈ L 2 . The grouping of solution components in a vector is underscored to distinguish it from an ordinary vector (coordinate representation).
2. Modeling concepts for structured systems. In this section, we briefly present some basic concepts for characterizing properties in a structured system. The concepts are applicable on the continuous level of our model problem as well as on the space-discrete level of its approximations. The core is the port-Hamiltonian framework that encodes geometric structures in an algebraic way. The geometric structures typically reflect fundamental physical properties, hence it is desirable to preserve them throughout all approximation steps. The partial Legendre transformation allows us to systematically investigate certain variable transformations related to the Hamiltonian of the system.

Port-Hamiltonian systems.
In recent years, the port-Hamiltonian framework has been found in many different areas of applications, see, e.g., [48,64,31,9,40]. Central is the Hamiltonian that is often related to an energy or entropy function in applications. In the finite-dimensional setting, H : R n → R is typically assumed to be convex and continuously differentiable. An important class of port-Hamiltonian systems with state z ∈ R n reads as follows, [64,9]: Here, the matrixJ(z) is anti-symmetric (i.e.,J(z) = −J(z) T ),R(z) is symmetric positive semi-definite (i.e.,R(z) =R(z) T has non-negative eigenvalues), and K ∈ R n×p . The choice of closing conditions depends on the application. Given an input u : [0, T ] → R p , the system (2.1) can, e.g., be closed by the equations e(t) = u(t), or f (t) = u(t). In the first case, the system reduces to an ordinary differential equation for z, whereas the second case yields a differential-algebraic equation of index 2.
In the port-Hamiltonian wording, z is called the energy variable, ∇H(z) the effort variable, and e and f the boundary effort and boundary flow, respectively. The system structure readily implies, among others, that the Hamiltonian is dissipated over time up to exchange with the boundary, i.e., for t ≥ 0 2.2. Partial Legendre transformation. The Legendre transformation plays an important role in classical mechanics, as it represents a link between the Lagrangian and the Hamiltonian modeling framework [68]. Further, there exists a rich theory on the implication of the Legendre transformation for duality principles, which are, e.g., used in nonlinear optimization [60]. In this paper we deal with a slight generalization of the standard approach, the so-called partial Legendre transformation [59]. Accordingly, we consider the following partitioning of a vector, z = [z 1 ; z 2 ] ∈ R n with z i ∈ R ni , n = n 1 + n 2 .
Definition 2.1 (Partial Legendre transformation). Let S = S 1 × S 2 ⊂ R n be a convex set with S i ⊂ R ni , n = n 1 + n 2 . The partial Legendre transformation of the function h : S → R with respect to the second (sub-vector) component is defined as For differentiable functions a coordinate transformation can be associated to the Legendre transformation.
By considering a 2 → g([z 1 ; a 2 ]) as a function parametrized inz 1 , the function h can be characterized as the partial Legendre transformation of g by standard results on Legendre transformations, cf. [59,60]. 3. Model problem. A class of nonlinear flows on networks is described by our model problem that covers, e.g., the barotropic Euler equations or electromagnetic waves. We write it in a specific form using a variable transformation that is related to a partial Legendre transformation. This is advantageous since parts of the coupling conditions become linear, which significantly simplifies the handling of network aspects. Moreover, it is key in the analysis of a variational principle. We show that the variational principle encodes an energy bound, local mass conservation and inherits a gradient structure strongly connected to the Hamiltonian of our model problem.
The nodes are grouped into interior nodes N 0 ⊂ N and boundary nodes N ∂ = N \N 0 , for an illustration of the network notation see Fig. 1.
On the network, function spaces are constructed by compositions of standard Sobolev spaces for every edge. The spatial domain is given as the union of edges Ω = {x : x ∈ ω, for ω ∈ E}. Note that every edge ω can be identified with an interval (0, l ω ) which is tacitly employed in the upcoming expressions. The space of squareintegrable functions on E reads L 2 (E) = b : Ω → R with b |ω ∈ L 2 (ω) for all ω ∈ E , where the subscript . |ω indicates the restriction of a function to the edge ω. The respective scalar product and norm read b,b = ω∈E ω bb dx and ||b|| = b, b for b,b ∈ L 2 (E). The weak (broken) derivative operator for functions on the network is defined by (∂ x b) |ω = ∂ x b |ω for ω ∈ E. The space of functions with square-integrable weak broken derivative is given as for all ω ∈ E denotes the space of piecewise smooth functions for k ≥ 0. The domains for boundary and coupling conditions are the sets of nodes N ∂ and N 0 . For b ∈ H 1 pw (E) we indicate node evaluations with squared brackets, i.e., b |ω [ν] ∈ R for ν ∈ N . They are well-defined by means of the trace theorem, [11]. Note that b |ω [ν] and b |ω [ν] may in general differ for ω =ω ∈ E. Following [24], a Right: Topology for gas pipeline network used in the numerical studies (Section 7.2). The larger red nodes ν i , i = 1, ..., 6, correspond to boundary nodes. The spatial domain of the pipes ω j , j = 1, ..., 8, i.e., the path from ν 1 to ν 3 , is highlighted by a dotted brown line.

Sobolev space incorporating a certain coupling condition at inner nodes is defined as
The boundary nodes N ∂ = {ν 1 , . . . , ν p } are assumed to be connected to exactly one edge each. Thus, a boundary operator T :

Strong form and variable transformation.
As model problem we consider a class of prototypical nonlinear partial differential equations for a flow on a network: The statez = [z 1 ; The solution components are interconnected by the coupling conditions for ω,ω ∈ E(ν), at ν ∈ N 0 . The system is closed by initial and boundary conditions of the formz for t ∈ [0, T ], ν ∈ N ∂ . We particularly assume one boundary data u ν : [0, T ] → R per boundary node ν ∈ N ∂ to be given, wheret ν denotes an appropriately chosen real-valued function. We refer to h as the Hamiltonian density and to as the (total) mass and the Hamiltonian, respectively. If (3.1) has a strong solution, the following two relations can be shown In our application (3.2a) has the interpretation of mass conservation, and (3.2b) is referred to as energy dissipation andz as energy variable from now on. Note that the structural properties also depend on the use of appropriate coupling conditions, in our case (3.1b) relates to conservation of the mass and the Hamiltonian at inner nodes, cf. [58,52].
Throughout our analysis, we assume that (3.1) has a smooth and unique solution. This assumption typically holds in the context of gas networks, which are governed by friction-dominated regimes [23,33,61], i.e., larger(·). Moreover, in many applications there is an additional dissipation term, see Appendix A. Our approximation scheme is constructed such that it mimics (3.2) on a discrete level, which enhances its stability in comparison to standard methods. The type of boundary conditions has no special influence on the proposed approximation ansatz, but many choices can be included very naturally.  and its inverse are continuously differentiable.
Proof. Theorem 2.2 is applicable, which implies that is the inverse ofẑ. By Assumption 3.1 the derivative d dzâ (z) is continuous, and a direct calculation shows that its eigenvalues λ 1 = 1 and λ 2 = ∇ 22 h(z) are strictly positive due to the positive definiteness of ∇ 2 h(z). Therefore, the derivative ofẑ is also a continuous function and reads , which follows from applying the implicit function theorem onâ locally for each z ∈ S. Because of Assumption 3.1, our model problem (3.1) is of strictly hyperbolic type. The underlying system matrix has one positive and one negative eigenvalue, such that one boundary condition per boundary node is required for a well-posed setup. We refer to [43,42] for details. Change of variables have played a crucial role in the theoretical and numerical analysis of hyperbolic systems, see e.g., [53,36,15,39]. The variable-transformed formulation considered in this work is related to the partial Legendre transformation of the Hamiltonian density.
(supplemented with closing conditions according to (3.1c)), where r(ā(t)) =r(ẑ(ā(t))) is non-negative and g the partial Legendre transformation of h. Thenẑ(ā) fulfills (3.1a)-(3.1b), and the Hamiltonian H of the system is given by Note that here and in the following, the solution is interpreted as a function in time with values in a function space. Moreover, by slight abuse of notation, the point- remark 3.4. The Hamiltonian H as function ofā is characterized in terms of the Legendre transformation g in (3.3c), cf. Lemma 2.3. The more direct characterization in the energy variable with respect to the Hamiltonian density h reads H(ā) =H(ẑ(ā)) = h(ẑ(ā), 1 . Note that Assumption 3.1 implies structural properties, e.g., strict convexity, onH rather than on H. remark 3.5. Due to the variable transformation the first coupling condition in (3.1b) becomes linear in the variable a 2 , although the model problem is of a general nonlinear form. This simplifies the handling of the network aspects, since we can use linear ansatz spaces in our approach, similar to [21]. Note that the incorporation, and beforehand, the identification of consistent coupling conditions for our model problem is a non-trivial issue, cf. [58,52].

Variational principle.
where the functions e, f : [0, T ] → R p with one entry for each boundary node ν i ∈ N ∂ , p = |N ∂ |, satisfy Proof. By testing the second equation of (3.3a) with b 2 ∈ H 1 div (E), integrating it over one edge ω = (ν,ν) ∈ E and using integration by parts, we obtain Here, the subscript . ω indicates the restriction onto the edge ω. Repeating the calculation for all edges and then summing up all equations, the interface terms at the inner nodes drop out, as b 2 ∈ H 1 div (E). This gives the second equation of the variational principle. The other equations follow similarly.
We use the variational principle to show local mass conservation and the energy dissipation equality.
Proof. The first relation, the local mass conservation, follows from testing the first equation of the variational principle (Theorem 3. The second equation, the global mass conservation, results from summing up the individual masses on the edges. As a preliminary step for the energy dissipation equality, we prove the existence of ξ(ā) for solutionsā of the variational principle with ξ(ā(t)) ∈ L 2 (E) such that Asā fulfills the variational principle, −∇ 1 g(ā), ∂ x b 2 has to be bounded for all b 2 ∈ H 1 div (E) and t ∈ [0, T ] due to the second equation of Theorem 3.6. The compatibility with ∂ x ζ = ξ follows from the definition of the broken derivative ∂ x . This shows, in turn, , which that b 1 → ∇ 1 h(ẑ(ā)), b 1 is an element of the dual space of L 2 (E). Thus, the existence of ξ(ā) fulfilling (3.4) follows from the Riesz representation theorem, as L 2 (E) with ·, · is a Hilbert space.
Let us now turn to the proof of the energy dissipation. A formal application of the chain rule leads to Using the variational principle withb = [ξ(ā); a 2 ] yields due to the non-negativity of r, which finishes the proof.
remark 3.8. The proof of the energy dissipation relies on the fact that for any functionā fulfilling the variational principle (Theorem 3.6), [−∇ 1 g(ā); a 2 ] can be identified with the gradient of the HamiltonianH(z) |z=ẑ(a) = h(z) |z=ẑ(a) , 1 with respect to the energy variablez. In the infinite-dimensional case, gradients of functionals, i.e., identifications of their derivative with an appropriately smooth function, do not exist for all choices of function spaces, cf. [68,12]. The gradient structure is a nontrivial property that our special variational principle has. It plays a fundamental role for the structure-preserving properties of our proposed approximation approach.
For each boundary node ν i ∈ N ∂ , we assume a boundary data u νi : [0, T ] → R to be given and pose a boundary condition of the form k i (e i , f i , u νi ) = 0. This equation only depends on the boundary effort e i and the boundary flow f i (Theorem 3.6). For a more concise notation, we collect all boundary conditions in a vector-valued function k : Note that the type of boundary condition affects the necessary regularity assumptions on the boundary data. To obtain for example a continuously differentiable solution a, it can be seen that a boundary condition of the type f i = u νi requires at least continuously differentiable data u νi , whereas for a boundary condition of the type e i = u νi continuity can be sufficient. In addition, the initial condition of the model problem are stated asā(0) =ā 0 .
4. Spatial approximation approach. The introduced variational principle is the basis for our structure-preserving spatial approximation approach which covers Galerkin projection and additional complexity reduction of the nonlinear terms. Local mass conservation and an energy bound for the approximations as well as a gradient structure related to the Hamiltonian of the system are kept under mild assumptions.
4.1. Galerkin approximation. The Galerkin approximation consists of a projection onto a suitable finite-dimensional ansatz space. It is applicable to a finite element discretization, but also in a more general setting, such as projection-based model reduction. Similar to the previous works [24,26,25] for linear and semilinear variants of our model problem, we impose few compatibility conditions as the only restriction.
In this paper only condition (A1) is explicitly used. Nonetheless, (A2) should be included as the analysis of the related steady state problem reveals, see [24,25] for the linear and semilinear case.
The proposed Galerkin approximation for our model problem reads as follows. System 4.3. Assumption 4.1 is supposed to hold. Given initial and boundary data,ā 0 ∈ V and u ν : for all b 1 ∈ V 1 , b 2 ∈ V 2 and fulfillā(0) =ā 0 and k(e, f , u) = 0.
The compatibility conditions on our ansatz spaces (Assumption 4.1) guarantee several structural properties for our Galerkin approximation.
Proof. Ifā fulfills (4.1), the expression −∇ 1 g(ā(t),b is well-defined for anyb ∈ ∂ x V 2 . Thus, it is well-defined for anyb ∈ V 1 due to the compatibility condition , b 1 is an element of the dual space of V 1 , which can be identified with a ξ(ā) ∈ V 1 by the Riesz representation theorem. The claim on the existence ofξ(ā) follows similarly by using that V 2 with ·, · is a Hilbert space.
Thus, our approximation inherits the gradient structure of the model problem. Here, in the finite-dimensional setting, the result is slightly stronger than in the infinitedimensional setting, where H 1 div (E) with ·, · is not a Hilbert space, cf. Remark 3.8.
Proof. The first relation implies that mass conservation holds in a pointwise sense, i.e., locally. To show it, note that by construction of the ansatz spaces ∂ t a 1 (t), Subtracting the right-hand side reveals that the L 2 -norm of ∂ t a 1 (t) + ∂ x a 2 (t) is zero, i.e., ∂ t a 1 (t) = −∂ x a 2 (t).
The proof of the energy dissipation equality is the same as that of Theorem 3.7, except for replacing the spaces L 2 (E) and H 1 div (E) by V 1 and V 2 and using Lemma 4.4. There are a few possibilities for constructing a suitable finite-dimensional approximation space V in the sense of Assumption 4.1, cf. [25,26]. We generate one by means of mixed finite elements. We divide the edges ω ∈ E into sub-parts T ω,k , k = 1, . . . , J ω and consider the polynomial spaces The sub-parts T ω,k induce a partitioning T E of the full network, on which we define the piecewise polynomial spaces Particularly, P q (T E ) inherits the coupling conditions encoded in H 1 div (E), which predefines one degree of freedom per inner node ν ∈ N 0 . The choice V 1 = Q q (T E ) and V 2 = P q+1 (T E ) for any q ≥ 0 yields a pair of compatible spaces.

Complexity reduction. System 4.3 features nonlinearities in
the balance equation (4.1b) due to the partial Legendre transformation g and the model function r. Depending on the model problem and the choice of the Galerkin ansatz space V, an additional complexity reduction can become necessary for an efficient numerical realization. We propose a quadrature-type approximation that preserves the structural properties. In particular, we allow for inexact integration, which is, e.g., crucial in the design of online-efficient reduced order models, cf. [27,4], but also may be convenient for high-order finite elements.
Given appropriate quadrature points x i ∈ Ω and quadrature weights w i for i ∈ I with I describing an index set, we introduce the complexity-reduced approximations ·, · c and || · || c of the L 2 -scalar product ·, · and L 2 -norm || · || as where b[x i ] refers to point evaluations of b. The following assumption ensures that the quadrature rule is well-defined and stable and that || · || c is a norm on V 1 and V 2 .
Assumption 4.6. For i ∈ I, the evaluations b[x i ] are well-defined for b ∈ V 1 ∪ V 2 , and the quadrature weights are positive, w i > 0. Moreover, there exists a constant C > 0 such that Our complexity reduction of System 4.3 modifies only the nonlinear balance equation (4.1b), it is replaced by Further, we consider the adjusted functional G c :ā → g(ā), 1 c associated to the partial Legendre transformation and the complexity-reduced Hamiltonian H c (ā) = ∇ 2 g(ā)a 2 − g(ā), 1 c . The gradient structure is kept, since the derivative of G c can be identified with [ξ(ā(t));ξ(ā(t))] ∈ V for t ∈ [0, T ] in the sense that ξ(ā(t)), b 1 = −∇ 1 g(ā(t)), b 1 c , ξ (ā(t)), b 2 = ∂ t ∇ 2 g(ā(t)), b 2 c holds for b 1 ∈ V 1 , b 2 ∈ V 2 . Note that b 1 → −∇ 1 g(ā), b 1 c and b 2 → ∂ t ∇ 2 g(ā), b 2 c are bounded linear functionals, i.e., elements of the dual space of V 1 and V 2 , respectively. This holds for any choice of inner product, in particular also for the L 2 -inner product used above. The energy dissipation equality becomes The proof of the structural properties follows straightforward the arguments of Lemma 4.4 and Theorem 4.5 using Assumption 4.6 on the equivalence of the norms. The local mass conservation is not affected by the complexity reduction at all. remark 4.7. Methods from literature on model reduction for port-Hamiltonian systems typically employ interpolation-type complexity reduction by DEIM. This yields either approximations, which are only port-Hamiltonian up to an approximation error [2,3], or the Hamiltonian structure is enforced by an additional symmetrization step [17], which sacrifices fidelity. For comparisons with our structure-preserving approach we refer to our follow-up paper [46].
holds by the chain rule. The matrix M is assumed to be symmetric positive definite. The weighted gradient and its sub-blocks relate to gradients with respect to the inner products induced by M and M i , respectively, cf. [1,49]. By choosing M as the mass matrix we can express the gradient structure observed in Lemma 4.4 in the algebraic setting.
Corollary 5.1. Let the following system matrices be given then System 4.3 can be equivalently stated as: Given initial and boundary data, a 0 ∈ R n and u ν : The expression relates to the term J T ∇ M 1 G(a(t)) in (5.2). Analogously, the structured algebraic representation of the integral expression with ∇ 2 g(ā) can be derived. Note that the block matrices due to the strict convexity of g with respect to its second component as consequence of Assumption 3.1. The bijectivity ofẑ follows then from Theorem 2.2. The port-Hamiltonian system is obtained from (5.2) via this coordinate transformation.

5.2.
Representation of complexity reduced system. An equivalence transformation of our approximation to a port-Hamiltonian system in standard form and a counterpart to Theorem 5.2 also exist for our complexity reduced system. Because of the strong similarity to the case without complexity reduction, we only highlight the differences. The structured coordinate representation belonging to the complexity reduction equals (5.2) except of the expressions G and R that are replaced by The respective Hamiltonian is H c (a) = ∇ 2 G c (a) · a 2 − G c (a). The coordinate transformation into standard port-Hamiltonian form is given byẑ c (a) = [a 1 ; ∇ M 2 G c (a)] for a = [a 1 ; a 2 ]. Its bijectivity can be concluded analogously as in the proof of Theorem 5.2. Note that for the strict convexity of Φ : a 2 → G c ([ā 1 ; a 2 ]) for fixedā 1 , Assumption 4.6 is needed. We can argue that forȃ 2 =ã 2 is valid, because the quadrature weights are positive (w i > 0) and there exists at least one quadrature point usingȃ 2 −ã 2 ∈ V 2 and the equivalence of || · || c and || · || on V 2 . The strict convexity of g yields (i.e., strict positivity at x j ) and ensures non-negativity at all other quadrature points x i , i ∈ I.
6. Time discretization. The underlying Hamiltonian structure of our space approximation can be exploited for analysis and control purposes but also for the derivation of energy-stable or energy-preserving time discretization schemes. Probably the simplest energy-stable scheme is provided by the implicit Euler-type method discussed in the following.
We consider an equidistant time grid t k = k∆ t ∈ [0, T ], k = 0, ..., K, K = T /∆ t with grid size ∆ t and indicate the temporal approximations by a respective superindex, e.g., a k ≈ a(t k ). System 6.1 (Implicit Euler-type scheme). Given initial data a 0 ∈ R n and boundary data u k νi ∈ R for ν i ∈ N ∂ , find a k = [a k 1 ; a k 2 ] ∈ R n and f k = [f k ν1 ; . . . ; f k νp ], e k = [e k ν1 ; . . . ; e k νp ] ∈ R p for k > 0 by solving with closing conditions a 0 = a 0 and k(e k , f k , u k ) = 0. System 6.1 can be interpreted as an implicit Euler discretization in the energy variable z and thus is of first-order convergence in time. A time-discrete counterpart to the energy bound from Theorem 4.5 can be derived. To do so, we make use of the next auxiliary result, which follows from the strict convexity of the Hamiltonian with respect to the energy variable.
Lemma 6.2. For the Hamiltonian H and the functional G given as in Corollary 5.1, it holds for a =ā that Proof. By Theorem 5.2, the coordinate transformationẑ : A n → R n to the energy variable is well-defined. The functionalH : R n → R,H(ẑ(a)) = H(a) thus represents the Hamiltonian as a function of the energy variable. As it is strictly convex by construction (due to Assumption 3.1), it holds The claimed inequality is equivalent to the latter one with z =ẑ(a),z =ẑ(ā), as can be seen by further employing that ∇H(z) |z=ẑ(a) = [−∇ 1 G(a); a 2 ], cf. Lemma 2.3.

Theorem 6.3 (Energy-dissipation inequality).
For any solution of System 6.1, it holds for 0 ≤ ℓ < k ≤ K and a ℓ = a k , Proof. Clearly, it is sufficient to consider the case ℓ = k − 1, as ℓ < k − 1 can be directly concluded from it. By Lemma 6.2 it holds Testing System 6.1 with b = [−∇ 1 G(a k ); a k 2 ] and using ∇ M i G(a k ) = M −1 i ∇ i G(a k ) for i = 1, 2, and K T [ã 1 ; a k 2 ] = K T 2 a k 2 = f k (which holds independently ofã 1 ), it follows due to the positive semi-definiteness of R(a). Setting together the two inequalities yields the assertion.
The proposed time discretization is energy-dissipative. In contrast to the energy dissipation equality in the time-continuous case (Theorem 4.5) we face here a dissipation inequality. The difference is due to the numerical dissipation of the implicit Euler-type method, which is D(t k ) = H(a 0 )−H(a k )+∆ t ( k j=1 e j ·f j −(a j 2 ) T R(a j )a j 2 ). Certainly, more sophisticated energy-stable or energy-preserving schemes can also be derived under the use of the underlying Hamiltonian structure, cf., [43,50,35,18], but this is beyond the scope of this paper.

7.
Application. The one-dimensional barotropic Euler equations and various simplifications of them are important representatives of our model problem. They are used, e.g., for gas network simulations, cf., [43,25,52]. In this section we demonstrate the feasibility of our approximation approach at the examples of an undamped dambreak test case with a shock and of a larger network in a friction-dominated regime relevant for pipelines in a gas transport network. For the last example, weighted edges must be considered which requires corresponding minor adjustments to the presented results, as we will briefly comment.

Barotropic Euler equations.
We consider the barotropic Euler equations with an additional friction term as an example for our model problem. Flow density and velocity ρ, v : [0, T ] × Ω → R are governed by where the pressure p is only a function of density. Furthermore, λ describes the friction factor that might be state-dependent and D the constant diameter of the pipe. In networks the circular cross-sectional pipe area A = π/4 D 2 usually acts as a weighting term for the respective edge ω ∈ E, see Appendix B. Following [55,66], we introduce P as a pressure potential that is characterized by the relation P ′′ (ρ) = p ′ (ρ)/ρ. Examples for compatible pairs of p and P are given in the test cases. Using the entropy as Hamiltonian H, system (7.1) can be formally rewritten in our standard form (3.1) for the energy variablez = [ρ; v], i.e., with Hamiltonian density h([ρ; v]) = ρv 2 /2+P (ρ) and non-negative friction-associated functionr([ρ; v]) = λ|v|/(2Dρ). The coupling conditions for ν ∈ N 0 take the form for ω,ω ∈ E(ν). The coupling conditions are entropy-preserving. Note that it is not trivial to find the entropy-preserving conditions without energy-based modeling or other sophisticated analytic tools, cf., [58,21,52]. The system is closed by appropriate initial and boundary conditions that are specified in the test cases. In our approach, the system is parametrized in the variables z 1 = ρ and ∇ 2 h(z) = ρv, i.e.,ā = [ρ; m] with m = ρv denoting the mass flux. In the hyperbolic theory these variables are a natural choice for the barotropic Euler equations that is typically used when a conservative formulation is discretized, cf., [42,66,20,58]. Although we do not approach the problem from the hyperbolic point of view, we use the variables as we rely on the port-Hamiltonian framework. Our proposed space discretization reads as follows. with r([ρ; m]) = λ|m|/(2Dρ 2 ) hold for all b 1 ∈ V 1 and b 2 ∈ V 2 . The system is closed with initial and boundary conditions,ā(0) =ā 0 and k(e(t), f (t), u(t)) = 0 for t ∈ [0, T ], for givenā 0 ∈ V and u ν : [0, T ] → R for ν ∈ N ∂ . In the case of complexity reduction, the inner products ·, · in the second equation are replaced by ·, · c . remark 7.2. Note that the preservation of port-Hamiltonian structure is a sufficient but not a necessary condition for a discrete energy bound (as in Theorem 4.5) to hold. A similar Galerkin approximation for the barotropic Euler equations that fulfills an energy bound but does not preserve this structure can be found in [21]. The underlying variational principle differs in the second equation, which reads in [21] In the reference no complexity reduction is considered and the derivation of an energystable time discretization seems to be more involved due to the loss of port-Hamiltonian structure. Note that port-Hamiltonian structure is also desirable for other reasons, e.g., for relative error estimates [30], or when using control theory [47,64,40].

Numerical studies.
Our numerical studies are performed with mixed finite elements in space and the proposed implicit Euler-type time discretization (System 6.1) with a uniform step size. The resulting nonlinear algebraic systems are solved by the Newton's method. When using finite elements of the lowest order, i.e., V 1 = Q 0 (T E ), V 2 = P 1 (T E ), the nonlinear integrals can be evaluated analytically, as they consist of polynomial expressions for the example. Higher-order finite elements are supplemented with complexity reduction by Gaussian quadrature of sufficiently high degree for Assumption 4.6 to hold. All results have been generated using MATLAB Version 9.1.0 (R2016b) on an Intel Core i5-7500 CPU with 16.0GB RAM. For better reproducibility, the code and the benchmark data are provided in [44]. 7.2.1. Dam-break test case. The undamped benchmark example on one edge is taken from [42,21] and referred to as dam-break problem. We use it here as a numerical stability test for our approach and therefore purposely omit to add a shock-capturing mechanism. For the latter we refer to, e.g., [63,42] and note that (numerical) dissipation can be included in our analysis without difficulty, see Appendix A. The model parameters are chosen as λ = 0, A = 1 and p(ρ) = 0.5ρ 2 , which corresponds to the isentropic Euler equations with pressure potential P (ρ) = p(ρ). The spatial domain is Ω = [−5, 5], and the initial and boundary conditions read ρ(0, x) = 2 − sgn(x), m(0, x) = 0 x ∈ Ω, and m(t, −5) = m(t, 5) = 0, t ≥ 0. As end time we take T = 2. The initial state has a discontinuity in the density ρ at x = 0, which propagates from left to right as a shock wave, whereas a rarefaction wave moves in the opposite direction, see Fig. 2 for an illustration.
For this setup, we compare our method with the finite element discretization proposed in [21] (cf., Remark 7.2), using the spaces V 1 = Q 0 (T E ) and V 2 = P 1 (T E ), see Fig. 3. Both space discretizations show a first-order convergence outside the region of the shock (x ∈ [−5, 0)) and overall very similar approximation quality. But our method is less dissipative, e.g., for ∆ x = 0.05 an energy loss of less than 1.4% is observed, i.e., H(ā(2)) ≈ 0.986 H(ā(0)), in contrast to the other one with a loss of more than 1.6%.
Further, we use the dam-break problem to showcase the entropy stability for higher-order discretization combined with inexact integration. Figure 4 exemplarily illustrates the results for the end time T = 2 and the choice V 1 = Q 3 (T E ), V 2 = P 4 (T E ) supplemented with a five-point Gaussian quadrature of the nonlinear integrals. The size of the elements ∆ x is varied, otherwise the same parameters as above are used. As to be expected due to the Gibbs phenomenon [63,6], oscillations occur near the shock when the solution is under-resolved, which is the case for ∆ x = 0.1 in this example. Nonetheless, the discrete solution stays bounded because of the entropy stability our method inherits. For the smaller grid sizes ∆ x = 0.02 and ∆ x = 0.01, the numerical dissipation coming from the time discretization seems to be sufficient to filter out the high frequency oscillations almost completely, as can be seen in the enlarged image , which implies the pressure potential P (ρ) = RT ρ log ((1 − RT αρ)ρ ⋆ /ρ). The friction factor λ is set constant over the whole network. With λ = 0.01 we choose a small, but typical value for gas networks. For illustration purposes we additionally consider λ = 0.003 being much smaller than observed in gas network modeling. The used pipeline network is visualized in Fig. 1. It consists of 38 pipes with diameters D ω between 0.4 and 1 [m] and lengths l ω between 5 to 74 [km]. The total pipe length of the network is 1008 [km]. The topology is a slight modification of [62,, whereby the compressors are replaced by pipes. Consequently, there is no compensation of the energy loss related to friction inside of the network. For the realization of network simulations with active elements, a coupling of the sub-components over their boundary efforts and flows can be used. We refer to [45] for details and to [44], where a simulation script for a network with a compressor is provided. The circular cross-sectional pipe areas A ω act as edge weights for ω ∈ E (cf. Appendix B). At the six boundary nodes ν i , i = 1, ..., 6 we prescribe the following boundary conditions ρ(t, ν 1 ) = (65 + u(t)) ρ ⋆ , ρ(t, ν 2 ) = (50 + u(t)) ρ ⋆ , ρ(t, ν 4 ) = (60 − u(t)) ρ ⋆ , ρ(t, ν 5 ) = 60 ρ ⋆ , ρ(t, ν 6 ) = 45 ρ ⋆ , Am(t, ν 3 ) = −100 (Am) ⋆ , with the time-varying input u, with reference time t ⋆ = 1 [h]. As initial condition the stationary solution belonging to the boundary conditions at t = 0 is taken. The end time is T = 5 t ⋆ . The space discretization is performed by means of the lowest-order finite elements with a uniform mesh on each pipe. The maximal length of a finite element satisfies ∆ x ≤ 200 [m], which yields 10 156 degrees of freedom in space. The time step is chosen as ∆ t = 1 [s], i.e., in total 18 000 steps for the simulation time of 5 hours. For this setup, the runtime on our machine is about 35 minutes per simulation. The computational cost is dominated by the nonlinear solves needed in each time step. The resulting density and mass flow are presented along the pipes ω j , j = 1, ..., 8, that build a path from the supplier node ν 1 to the consumer node ν 3 (see Fig. 1). To illustrate the impact of the damping on the temporal evolution, the solution at fixed space points along the path is plotted in Fig. 5. Stronger damping effects and lower peaks are observed for the larger, realistic friction factor (λ = 0.01). Moreover, the sharp input profile u is more clearly transferred to the mass flow, whereas the dynamics of the density is smoothed out stronger. The simulation time the moving profile needs to travel from the end of pipe ω 2 to the end of pipe ω 7 is about 0.5t ⋆ for density and mass flow. For t ≥ 3t ⋆ the state seems to be almost at rest in the considered spatial domain for both choices of λ, whereby the equilibrium is reached slightly faster for the model with the application-relevant friction factor (λ = 0.01). For this model, we also visualize the solution in space in Fig. 6, using time snapshots along the path. The pipe junctions are indicated by vertical lines. Note that the mass flow is discontinuous at junctions with more than two pipes, which is in accordance to the coupling conditions (3.1b). As the spatial representation illustrates, no shocks are to be expected in the considered friction-dominated regime, and the solution has a rather simple structure in space. The latter motivates the use of model reduction techniques to speed up gas network simulations, cf., [43,32,19], which in contrast to finite element methods rely on ansatz functions with global support over the full network.
The results of this paper lay the theoretical foundation for generalizing our model order reduction approach [26] for the linear damped wave equation to a model orderand complexity-reduction approach for a general nonlinear flow problem class. In our follow-up paper [46], the algorithmic aspects related to snapshot-based model reduction in our structure-preserving framework are investigated.
8. Conclusion and outlook. For many practically relevant dynamical systems, an underlying Hamiltonian structure can be identified. In this paper we proposed a structure-preserving approximation approach relying on a port-Hamiltonian formulation for a class of nonlinear flows on networks. Its basis is a variational principle that inherits the Hamiltonian structure. Further, the parametrization of the solution takes a prominent role. We take a mix of energy and co-energy variables and thus avoid to fully change to the co-energy variables (which are equal to the entropy variables for the Euler equations). To treat the parametrization and the nonlinearities in a systematic manner, we employ the theory on Legendre transformations. Apart from the port-Hamiltonian structure, our approximations ensure local mass conservation and an energy bound under mild assumptions. In particular, quite general Galerkin projections and complexity reduction of nonlinearities by a quadrature-type ansatz are covered by our analysis. We showcased the applicability and good stability properties of our approximation at the example of the barotropic Euler equations. However, our approach can also be applied to any p-system and other symmetrizable hyperbolic systems in two variables, e.g., in the context of electromagnetic waves.
The results of this paper can also be used as a theoretical basis for structurepreserving model reduction methods, e.g., for gas network systems. Note that the realization of compatible reduced models involves some non-trivial algorithmic issues. In parts, they are similar to the ones occurring in symplectic model order reduction [57,2,3], some are covered in our former works [43,25,26]. A thorough discussion of the algorithmic aspects related to snapshot-based model order-and complexityreduction in our framework can be found in our follow-up paper [46]. Another interesting direction for future research is the extension and adaption of our approach to other classes of model problems and more involved applications, e.g., to multidimensional systems with a similar underlying Hamiltonian structure [30,5,61], or thermal Euler equations with additional dissipation terms [54,37].
Data Availability. The MATLAB code used to generate the presented numerical results and some supplementary tests, including a simulation of a network with a compressor, can be found under the DOI 10.5281/zenodo.6372667, see [44].
Appendix A. Inclusion of dissipation.
Our model problem (3.1) is of hyperbolic type. This changes when dissipation effects in form of a second-order derivative in space are included. We assume them to be described by a non-negative term d : R 2 → R + , similar as in [21]. In generalization to (3.1), let the statē z : [0, T ] × Ω → R 2 be governed by z(t, x))∂ x ∇ 2 h(z(t, x))) with coupling conditions at ν ∈ N 0 given by ω∈E(ν) n ω [ν]∇ z2 h(z |ω (t, ν)) = 0, s(z |ω (t, ν)) = s(z |ω (t, ν)) for ω,ω ∈ E(ν), where s(z) = ∇ 1 h(z) − d(z)∂ x ∇ 2 h(z) forz = [z 1 ; z 2 ] ∈ C 1 pw (E). To close the system, initial conditions and one boundary condition per boundary node ν ∈ N ∂ have to be prescribed, similarly to the case without dissipation. Sufficiently smooth solutions can be shown to fulfill the energy dissipation equality Our spatial approximation approach is transferable to this problem. Respective discrete energy bounds and port-Hamiltonian structure can be shown with very minor adjustments to our derivations from the main part.
Appendix B. Edge weights. Pipelines in gas transport networks are typically modeled with cross-sectionally averaged dynamics. Hence, we include the cross-sectional pipe area A ω for ω ∈ E as edge weight in our approach. Note that the edge weighting affects all integral-related expressions and definitions from Section 3.1. The inner product ·, · becomes b,b = ω∈E A ω ω b[x]b[x]dx and the incidence mapping n ω [ν] = A ω for ω = (ν,ν) for someν ∈ N −A ω for ω = (ν, ν) for someν ∈ N , which modifies boundary and coupling conditions, and thus the boundary operator T : H 1 pw (E) → R p and the function space H div (E). Moreover, the Hamiltonian is altered to H(z) = h(z), 1 = ω∈E A ω ω h(z)dx.