Numerical approximation of port-Hamiltonian systems for hyperbolic or parabolic PDEs with boundary control
Andrea Brugnoli, Ghislain Haine, Anass Serhani, Xavier Vasseur
NNumerical approximation of port-Hamiltonian systemsfor hyperbolic or parabolic PDEs with boundarycontrol ∗ Brugnoli , Andrea † Haine , Ghislain ‡ Serhani , Anass § Vasseur , Xavier ¶ ISAE-SUPAERO, Université de Toulouse, France.10 Avenue Edouard Belin, BP-54032, 31055 Toulouse Cedex 4.
Abstract
We consider the design of structure-preserving discretization methods for the solution of systemsof boundary controlled Partial Differential Equations (PDEs) thanks to the port-Hamiltonian for-malism. We first provide a novel general structure of infinite-dimensional port-Hamiltonian systems(pHs) for which the Partitioned Finite Element Method (PFEM) straightforwardly applies. Theproposed strategy is applied to abstract multidimensional linear hyperbolic and parabolic systemsof PDEs. Then we show that instructional model problems based on the wave equation, Mindlinequation and heat equation fit within this unified framework. Secondly we introduce the ongoingproject
SCRIMP (Simulation and ContRol of Interactions in Multi-Physics) developed for the nu-merical simulation of infinite-dimensional pHs.
SCRIMP notably relies on the
FEniCS open-sourcecomputing platform for the finite element spatial discretization. Finally, we illustrate how to solvethe considered model problems within this framework by carefully explaining the methodology. Asadditional support, companion interactive
Jupyter notebooks are available.
Keywords : port-Hamiltonian systems; Partial differential equations; Boundary control; Structure-preserving discretization; Finite Element Method
Mathematical classification (2010) : 65M60; 35L90; 35K90
The efficient numerical simulation of complex multiphysics systems is ubiquitous in Computational Sci-ence and Engineering. Although a wide range of methods exists to tackle specific problems, they oftenlack of versatility and adaptability, especially when the modelling is of increasing complexity as in real-world applications.Infinite-dimensional port-Hamiltonian systems (pHs) have been first introduced in [55] using thelanguage of differential geometry. They provide a powerful tool to model complex multiphysics opensystems (whether or not being linear) for control purpose. A wide range of physical systems has beenwritten within this formalism, see e.g. [36, 57, 58]. This twenty-year-old framework [43] enjoys niceproperties, such as the relevant physical meaning of the variables, a useful underlying linear structure(namely Stokes-Dirac structure) which encodes the power balance satisfied by the Hamiltonian (oftenchosen as an energy), and last but not least: the interconnection of multiple pHs remains a pHs. Thisallows a “modular” modelling of complex multiphysics systems.Since then, many researchers have developed numerical methods to discretize these systems in astructure-preserving manner, hence keeping the advantages of the infinite-dimensional pHs. Such meth-ods aim at constructing an approximate finite-dimensional pHs at the discrete level. Our aim is toshow the versatility of PFEM thanks to a new unified framework, and to introduce the ongoing project ∗ This work is supported by the project ANR-16-CE92-0028, entitled
Interconnected Infinite-Dimensional systems forHeterogeneous Media , INFIDHEM, funded by the French National Research Agency (ANR) and the Deutsche Forschungs-gemeinschaft (DFG). Further information is available at https://websites.isae-supaero.fr/infidhem/the-project. † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] a r X i v : . [ m a t h . NA ] S e p CRIMP together with companion
Jupyter notebooks [17]. Each particular example discussed herehas been treated previously [14, 19, 46, 47]. Here, the existence of an underlying common structurefor many pHs is highlighted. Obtaining such a general scheme for infinite-dimensional pHs is of majorimportance for control purposes [51], and for coupling atomic elements into a more complex system withthe guarantee of well-preserved energy exchanges between subsystems [34].The first proposed structure-preserving scheme for pHs dates back to [25], where the authors proposeda mixed finite element spatial discretization for hyperbolic systems of conservation laws. Pseudo-spectralmethods relying on higher-order global polynomial approximations were studied in [41]. Unfortunatelythis method seems to be limited to the one-dimensional case. A finite difference method with staggeredgrids was developed in [52] for two-dimensional domains, but complex geometries are then difficult totackle. Weak formulations leading to Galerkin numerical approximations began to be explored in thepast few years. In [33] the prototypical example of hyperbolic systems of two conservation laws has beendiscretized by a weak formulation. However the construction of the necessary power-preserving mappingsis not straightforward on arbitrary meshes. All these methods require ad hoc implementations, and areusually restricted to particular cases of pHs. Furthermore, since they do not rely on well-establishedand versatile numerical libraries, using such techniques remains confined within a small community ofexperts. We refer the reader for a more complete overview of structure-preserving discretization for pHsto [43, 31] and the references therein.Thanks to [19] it has become clear that there exists a deep relation between structure-preservingdiscretization of pHs and Mixed Finite Element Method (MFEM). Indeed velocity-stress formulationsfor the wave dynamics [29] and elastodynamics problems are of Hamiltonian type and their mixeddiscretization preserves this structure for closed systems. This leads to the intuition that a MFEM maybe used to discretize the underlying geometric structure of pHs in a unified way, even for open systems,translating the infinite dimensional Stokes-Dirac structure into a finite Dirac structure. The discretizationstrategy relies on the partitioned structure of the problem and for this reason goes under the name ofPartitioned Finite Element method (PFEM). This method proves nice convergence properties, see e.g. [26] for a recent proof on the wave equation, that does not require the fulfilment of the usual inf–sup condition for MFEM, generalizing the results cited in [28, Remark 6].It has to be pointed out that the core idea of PFEM, i.e. performing an integration by parts ona partition of the weak formulation of the system of equations, has already been proposed for closed hyperbolic systems in [28]. Therein, the formulations are called either primal–dual or dual–primal ,depending on the chosen partition of the system.The major difference between MFEM and PFEM relies on the choice of the test functions in the weakformulation, hence on the finite element form functions. Indeed, in PFEM, they never carry homogeneousconditions. In e.g. [12, Section 7.1], it is shown that for a Dirichlet control, test functions are takenin the kernel of the Dirichlet trace. As already mentioned, in [28], the proposed primal–dual and dual–primal discretizations are then suitable for the structure-preserving discretization of closed systems.Nevertheless, by keeping these homogeneous conditions in the test functions, not only the application ofDirichlet control is difficult, but the definition of the Neumann observation, necessary for the discretepower balance, would be more complex. PFEM aims at easing the mimicking of the continuous powerbalance at the discrete level, by relaxing the test functions used in [28]. To the best of our knowledge,this relaxation has not been investigated yet in the case of boundary controlled wave-like systems, andthis probably comes from the fact that MFEM have been first developed for elliptic systems. Indeed, inthis case, such a boundary condition is mandatory for well-posedness (especially to obtain the ellipticityin the kernel condition [12, Eq. (5.1.7)]), while PFEM is made for evolution systems, and more especiallyfor pHs. Furthermore, we are not aware of applications of the studied scheme, being called MFEM orPFEM, to the boundary controlled heat equation, beside the first attempt presented in [46, 47].In our opinion, the driving forces of PFEM are threefold: first, PFEM takes collocated boundarycontrols and observations into account in a simple manner; secondly, PFEM is structure-preserving,meaning in particular that the discrete power balance perfectly mimics the continuous one; thirdly, theimplementation of PFEM only relies on existing finite element libraries, such as FEniCS [2], selectedin the ongoing project SCRIMP for its robustness and efficiency. Last but not least, the pHs pointof view allows us to separate axioms of physics (such as conservation laws) from constitutive laws andequations of state (such as Hooke’s law and ideal gas law). PFEM is based on this separation, providingthe possibility to tackle parabolic or nonlinear systems, at the price of solving a finite-dimensional port-Hamiltonian Differential Algebraic Equation (pHDAE) . In the particular case of linear hyperbolic PDE,as shown in Section 2, the constitutive laws can be easily ( i.e. without matrix inversions) taken intoaccount in order to recover an
Ordinary Differential Equation (ODE) .2FEM could also be named e.g. extended MFEM or relaxed MFEM . Since only evolution systems areconsidered (not necessarily of hyperbolic type, see e.g. Section 3.4), relaxed conditions for the selectionof the test functions hold, hence for the finite elements as well. We choose to follow the terminologyintroduced in [19] and widely used since then. Furthermore, it emphasizes the pH formulation of theinitial system to discretize.
Main contributions
We first aim at presenting the strategy of the structure-preserving discretizationPFEM, in a new unified abstract framework, allowing for an easy application to a wide class of boundarycontrolled partial differential equations. Then, in order to show the versatility of our approach, we suc-cessively apply PFEM to the boundary controlled wave equation, the boundary controlled Mindlin platemodel, and the boundary controlled heat equation with a thermodynamically well-founded Hamiltonian(namely the internal energy , instead of the quadratic functional commonly used). Taking advantage of thestrong underlying structure, we finally describe a unified object-oriented implementation of these models via
PFEM. Companion interactive
Jupyter notebooks [17] are discussed to illustrate our methodology.
Structure of the manuscript
The manuscript is organized as follows. In Section 2 the abstractpHs framework is introduced, with a particular focus on both hyperbolic and parabolic linear systemsof partial differential equations. In Section 3 the general structure-preserving discretization is presented,and then specialized on the two cases previously mentioned. In Section 4 the ongoing environmentSCRIMP is described in detail. In Section 5 the three companion interactive
Jupyter notebooks [17]are thoroughly explained. Conclusions and perspectives are finally drawn in Section 6.
In this section, we introduce an abstract class of pHs and their underlying geometric structure: the Diracstructure for the finite dimensional case and the Stokes-Dirac structure for the infinite dimensional case.For the infinite-dimensional case it is shown how hyperbolic and parabolic systems easily fit into thisframework.
State representation
Let us begin with a classical definition of a pHs in finite dimension. Considerthe time-invariant dynamical system [54]: d x dt = ( J ( x ) − R ( x )) ∇ H ( x ) + B ( x ) u , y = B ( x ) (cid:62) ∇ H ( x ) , (1)where H ( x ) : X (cid:39) R n → R , the Hamiltonian, is a real-valued function of the vector of energy variables x , bounded from below. Matrix-valued functions J ( x ) (the structure operator ) and R ( x ) (the dissipative or resistive operator ) are skew-symmetric and symmetric positive semi-definite respectively. The control u ∈ R m is applied thanks to the matrix-valued control function B ( x ) of size m × n . Variable y ∈ R m isthe power conjugated output to the input.Such a system is called a port-Hamiltonian system , as it arises from the Hamiltonian modelling of aphysical system and it interacts with the environment via the input u and the output y , included in theformulation. The vector ∇ H ( x ) is made of the co-energy variables .Due to the structural properties of J ( x ) and R ( x ) , the port-Hamiltonian system enjoys the nicefollowing power balance : dHdt = − ( ∇ H ( x )) (cid:62) R ( x ) ∇ H ( x ) + u (cid:62) y ≤ u (cid:62) y , (2)meaning that R ( x ) accounts for dissipation, and that the input–output product corresponds to the powersupplied to (or took from) the system, through the control u . Flow–effort representation
Consider two finite dimensional vector spaces E = F (cid:39) R n . The ele-ments of F are called flows , while the elements of E are called efforts . Those are port variables and theircombination gives the power flowing inside the system. The space B = F × E is called the bond spaceof power variables. Therefore, identifying E as the dual of F , the power is defined as (cid:104) e , f (cid:105) := e ( f ) .3 efinition 1 ([23], Def. 1.1.1) Given the finite-dimensional space F and its dual E with respect tothe inner product (cid:104)· , ·(cid:105) F : F × F → R , consider the symmetric bilinear form: (cid:104)(cid:104) ( f , e ) , ( f , e ) (cid:105)(cid:105) := (cid:104) e , f (cid:105) + (cid:104) e , f (cid:105) , where ( f i , e i ) ∈ B, i = 1 , . A Dirac structure on B := F × E is a subspace D ⊂ B , which is maximally isotropic under (cid:104)(cid:104)· , ·(cid:105)(cid:105) .Equivalently, a Dirac structure on B := F × E is a subspace D ⊂ B which equals its orthogonal companionwith respect to (cid:104)(cid:104)· , ·(cid:105)(cid:105) , i.e. D = D [ ⊥ ] , where: D [ ⊥ ] := (cid:110) ( f , e ) ∈ B | (cid:68)(cid:68) ( f , e ) , ( (cid:101) f , (cid:101) e ) (cid:69)(cid:69) = 0 , ∀ ( (cid:101) f , (cid:101) e ) ∈ D (cid:111) . The connection between the concept of Dirac structure and pHs in its canonical form (1) is achievedby considering the following ports :• the storage ports ( f x , e x ) := (cid:0) d x dt , ∇ H ( x ) (cid:1) ∈ R n × R n , made of the storage flow f x (time-derivativeof the energy variables) and storage effort e x (the co-energy variables);• the resistive (or dissipative) ports ( f R , e R ) ∈ R k × R k , made of the resistive (or dissipative) flow f R and resisitive (or dissipative) effort e R ;• the interconnection ports ( f u , e u ) := ( − y , u ) ∈ R m × R m , made of the interconnection flow f u and interconnection effort e u .Assuming that the matrix R ( x ) has constant rank, from classical matrix factorizations there existmatrices G (not necessarily square, of size k × n ) and K symmetric positive semi-definite of size k × k such that R = G (cid:62) KG . These notations at hand, the pHs (1) rewrites: f x ( x ) f R ( x ) f u ( x ) = J ( x ) − G ( x ) (cid:62) B ( x ) G ( x ) − B ( x ) (cid:62) e x ( x ) e R ( x ) e u ( x ) , (3)together with the (dissipative or resistive) constitutive relation : e R ( x ) = K ( x ) f R ( x ) . (4)It is clear that the extended structure operator J e appearing in (3) is skew-symmetric of size ( n + k + m ) × ( n + k + m ) . Its graph is a Dirac structure with respect to the Euclidean inner product, as a kernelrepresentation, see [54]. Hence, it comes: (cid:104) e x ( x ) , f x ( x ) (cid:105) R n + (cid:104) e R ( x ) , f R ( x ) (cid:105) R k + (cid:104) e u ( x ) , f u ( x ) (cid:105) R m = 0 . Noting that dHdt = (cid:104) e x ( x ) , f x ( x ) (cid:105) R n (by definition of the storage port) leads to: dHdt = − (cid:104) e R ( x ) , f R ( x ) (cid:105) R k − (cid:104) e u ( x ) , f u ( x ) (cid:105) R m , and thanks to the symmetry of R , (4) gives: dHdt = − (cid:104) f R ( x ) , K ( x ) f R ( x ) (cid:105) R m − (cid:104) e u ( x ) , f u ( x ) (cid:105) R m . Finally, from (3) and the definition of the storage and interconnection ports, the power balance (2) isrecovered.The relation between (1) and (3)–(4) can be understood as follows: the power balance (2) is encoded in the Dirac structure (obtained from the extended structure operator J e ) together with the resistiveconstitutive relation. Remark 1
The canonical Euclidean inner product has been used here, but other inner products are al-lowed to take into account mass matrices (symmetric positive definite) on the left-hand side of (1) , (3) ,and (4) . This is crucial after the spatial discretization procedure. This corresponds to a kernel represen-tation of Dirac structure [54]. System 1 is a pHs in canonical form. Recently, finite-dimensional differential algebraic port-Hamiltonian systems (pHDAE) have been introduced both for linear [9] and nonlinear systems [40].This enriched description shares not only all the crucial features of ordinary pHs, but also easily ac-counts for algebraic constraints, time-dependent transformations and explicit dependence on time in theHamiltonian. The application of the proposed discretization method naturally leads to pHDAEs. Indeed,a constitutive relation between f x := d x dt and e x := ∇ H ( x ) is needed to be well-defined. But PFEM takesinto account constitutive relations apart from (3) as constraints. However, as shown later in Sections 3.2and 3.3, the method simplifies in the case, for instance, of linear hyperbolic systems.4 .2 Infinite-dimensional port-Hamiltonian systems In this section an infinite-dimensional generalization of pHs is presented. For sake of readability, the(Stokes-)Dirac structure is first defined, and secondly, infinite-dimensional pHs are then described inboth hyperbolic and parabolic cases. A more general framework can be designed, but this goes beyondthe aim of this present work.
Structure operator
As to avoid functional difficulties, the analogue of the extended structure operatorwill not be written as in (3). More precisely, the control operator will not be included in an extendedstructure unbounded operator, but given apart. The Stokes-Dirac will be then obtained thanks to astructure operator related to the boundary control operator through an abstract Green formula. However,like in finite dimension, the aim is to establish a link between flow and effort variables. Most importantly,the underlying Stokes-Dirac structure must encode the power balance of the dynamical system understudy.Consider a Lipschitz domain Ω ⊂ R d , d ∈ { , , } , and the relation: (cid:18) f f (cid:19) = (cid:20) −L ∗ L (cid:21) (cid:18) e e (cid:19) , e ∈ H L : { e ∈ L (Ω , A ) | L e ∈ L (Ω , B ) } , e ∈ H −L ∗ : { e ∈ L (Ω , B ) | − L ∗ e ∈ L (Ω , A ) } . (5)By L (Ω , X ) we denote the space of square integrable X -valued functions. Symbol A , B denote either thespace of scalars R , vectors R d , symmetric tensors R d × d sym =: S or a Cartesian product of those, dependingon the particular example. The operator L : H L → L (Ω , B ) is a generic differential, and therefore linearbut unbounded, operator. The notation L ∗ : H −L ∗ → L (Ω , A ) denotes the formal adjoint of L , definedby the relation: (cid:104)L e , e (cid:105) L (Ω , B ) = (cid:104) e , L ∗ e (cid:105) L (Ω , A ) , e ∈ C ∞ (Ω , A ) , e ∈ C ∞ (Ω , B ) . (6)Of course, for (5) to be well-defined, constitutive relations are needed. Only physical laws will betaken into account when constructing the above relation on concrete examples. As pointed out in theintroduction, PFEM aims at both preserving this relation and discretizing constitutive laws to close thesystem. Remark 2
One can be confused by the lack of evolution in time in (5) . However, this emphasizes animportant paradigm in the proposed point of view: this relation translates the time-independent geometricstructure of the pHs as an equation with differential operator (ill-posed on its own), while constitutiverelations will bring back the time dependency of the problem. In particular, some flows must be the timederivative of the energy variables.
Throughout the paper, (cid:104)· , ·(cid:105) X denotes the inner product of the Hilbert space X . Definition (6) isanalogous to Definition 5.80 in [44]. In Section 3.2, the operator L is the gradient, denoted by grad ,and its formal adjoint is the divergence, denoted by div , from the so-called Green’s formula (integrationby parts). In Section 3.3, the operator L contains both grad and Grad . This latter corresponds to thesymmetric part of the gradient and represents the deformation tensor in continuum mechanics:
Grad( e ) := 12 (cid:0) ∇ e + ∇ (cid:62) e (cid:1) , e ∈ C ∞ ( R d ) . The formal adjoint of
Grad is the tensor divergence
Div . For a tensor field E : Ω → M := R d × d , withcomponents e ij , the divergence is a vector, defined columnwise as: Div( E ) := (cid:32) d (cid:88) i =1 ∂ x i e ij (cid:33) j =1 ,...,d . Finally, in Section 3.4, L is made of the gradient and the identity operator. Stokes-Dirac structure
Definition 1 still remains valid in infinite dimension. Nevertheless, as statedabove, the structure operator in (5) is not extended to include the control operator. Hence an additionalassumption has to be made for (cid:20) −L ∗ L (cid:21) to define a Dirac structure in relation with a pHs coming5rom boundary control of partial differential equations. In other words, a Stokes-Dirac structure requiresthe specification of boundary variables in order to express a general power conservation property for open physical systems. This assumption is based on the so-called Stokes’ theorem (also known as thedivergence theorem, Gauss’s theorem or Ostrogradsky’s theorem) and its corollaries, as the Green’sformula. Assumption 1 (Abstract Green’s formula)
The operator L is assumed to satisfy the abstractGreen’s formula: (cid:104)L e , e (cid:105) L (Ω , B ) − (cid:104) e , L ∗ e (cid:105) L (Ω , A ) = (cid:104) Γ e , Γ ⊥ e (cid:105) V ∂ , ( V ∂ ) (cid:48) , ∀ e ∈ H L , e ∈ H −L ∗ , (7) where the right-hand side is the duality bracket at the boundary, on a well-suited boundary functionalspace V ∂ for some trace operators Γ , Γ ⊥ . From now on, this duality bracket will be denoted by (cid:104)· , ·(cid:105) ∂ Ω with a slight abuse of notation. Remark 3
This abstract formula is well-known in the boundary control systems theory, see e.g. [53,Chapter 10].
Remark 4
In practice, equation (7) dictates the causalities , i.e. the possible choices for the boundarycontrol u ∂ and the boundary observation y ∂ , via the equality (cid:104) Γ e , Γ ⊥ e (cid:105) ∂ Ω = (cid:104) u ∂ , y ∂ (cid:105) ∂ Ω (with aslight abuse of notation for the right-hand side to make sense). Of course, the admissible causalities arealso related to the well-posedness of the system under study, and in particular to the definitions of theboundary functional spaces. For sake of simplicity, a focus on the two following causalities will be made. Let the boundaryvariables associated to system (5) be defined by: e ∂ = Γ ⊥ e ∈ ( V ∂ ) (cid:48) , f ∂ = − Γ e ∈ V ∂ , (8)or the other way: e ∂ = Γ e ∈ V ∂ , f ∂ = − Γ ⊥ e ∈ ( V ∂ ) (cid:48) . (9)In light of (7), systems: (cid:18) f f (cid:19) = (cid:20) −L ∗ L (cid:21) (cid:18) e e (cid:19) , (cid:18) e ∂ f ∂ (cid:19) = (cid:20) Γ ⊥ − Γ (cid:21) (cid:18) e e (cid:19) , (10)and: (cid:18) f f (cid:19) = (cid:20) −L ∗ L (cid:21) (cid:18) e e (cid:19) , (cid:18) e ∂ f ∂ (cid:19) = (cid:20) Γ − Γ ⊥ (cid:21) (cid:18) e e (cid:19) , (11)define Stokes-Dirac structures with respect to the bilinear pairing: (cid:10)(cid:10) ( f , f , f ∂ , e , e , e ∂ ) , ( f , f , f ∂ , e , e , e ∂ ) (cid:11)(cid:11) = (cid:10) f , e (cid:11) L (Ω , A ) + (cid:10) f , e (cid:11) L (Ω , B ) + (cid:10) f , e (cid:11) L (Ω , A ) + (cid:10) f , e (cid:11) L (Ω , B ) + (cid:10) f ∂ , e ∂ (cid:11) ∂ Ω + (cid:10) f ∂ , e ∂ (cid:11) ∂ Ω . Obviously, for systems (10) and (11) to be well-defined, constitutive relations are needed.In the remaining part of this section, only (10) will be treated in details. The other canonical causality(11), figuring in Section 3.4, straightforwardly follows with the same strategy.
Hyperbolic systems
In the hyperbolic case, both flows represent the dynamics of the independentenergy variables α , α . The Hamiltonian is a generic functional of these variables H = H ( α , α ) .The co-energy variables are by definition the variational derivatives (see e.g. [56]) of H with respect tothe energy variables: f = ∂ t α , e := δ α H, f = ∂ t α , e := δ α H. (12)Then system (10) possesses the equivalent state representation: (cid:18) ∂ t α ∂ t α (cid:19) = (cid:20) −L ∗ L (cid:21) (cid:18) δ α Hδ α H (cid:19) , (cid:18) u ∂ y ∂ (cid:19) = (cid:20) Γ ⊥ Γ (cid:21) (cid:18) δ α Hδ α H (cid:19) . (13)It holds u ∂ = e ∂ , y ∂ = − f ∂ . The power balance is naturally embedded in the Stokes-Dirac structuredefined by (10): ddt H ( α , α ) = (cid:104) y ∂ , u ∂ (cid:105) ∂ Ω . (14)6 inear hyperbolic systems The system is linear when the Hamiltonian has the form: H ( α , α ) = 12 (cid:104) α , Q α (cid:105) L (Ω , A ) + 12 (cid:104) α , Q α (cid:105) L (Ω , B ) , where Q , Q are positive symmetric operators, bounded from below and above: m I A ≤ Q ≤ M I A , m I B ≤ Q ≤ M I B , m > , m > , M > , M > , with I A and I B the identity operators in L (Ω , A ) and L (Ω , B ) respectively. In this case, the co-energyvariables are given by: e := δ α H = Q α , e := δ α H = Q α . (15)Since Q , Q are positive and bounded from below and above, it is possible to invert them to obtain: α = Q − e =: M e , α = Q − e =: M e , (16)giving rise to the co-energy formulation . The Hamiltonian is rewritten as: H ( e , e ) = 12 (cid:104) e , M e (cid:105) L (Ω , A ) + 12 (cid:104) e , M e (cid:105) L (Ω , B ) , (17)and a linear hyperbolic pHs (10) can be expressed as: (cid:20) M M (cid:21) (cid:18) ∂ t e ∂ t e (cid:19) = (cid:20) −L ∗ L (cid:21) (cid:18) e e (cid:19) , (cid:18) u ∂ y ∂ (cid:19) = (cid:20) Γ ⊥ Γ (cid:21) (cid:18) e e (cid:19) . (18)In this particular case, the constitutive relations needed for system (10) to be well-defined are givenby (15), and then directly included in (18). In Sections 3.2 and 3.3, it will be shown that PFEM leadsdirectly to a finite-dimensional pHs of the form (1) with R ( x ) = 0 . This simplification considerablyfacilitates the solution in time, as (1) is an Ordinary Differential Equation (ODE). Parabolic systems
In this case, the first flow f still represents a dynamics ∂ t α of the energy variable α . The Hamiltonian then reads H = H ( α ) , and its variational derivative gives the co-energy variable e = δ α H .The second flow f represents an extra flow related to the effort variable e appearing in the dynamicsof the energy variable α . The relation is given implicitly by a mapping G as G ( f , e ) = 0 . Then,pHs (10) of parabolic type is expressed as: (cid:18) ∂ t α f (cid:19) = (cid:20) −L ∗ L (cid:21) (cid:18) δ α H e (cid:19) , (cid:18) u ∂ y ∂ (cid:19) = (cid:20) Γ ⊥ Γ (cid:21) (cid:18) δ α H e (cid:19) , G ( f , e ) = 0 . (19)In Section 3.4, an example of a parabolic-type pHs (11) is studied. It will be shown that the PFEMstructure-preserving discretization of such a system naturally leads to a finite-dimensional pHDAE.Again, the power balance is naturally embedded in the Stokes-Dirac structure defined by (10): ddt H ( α ) = − (cid:104) f , e (cid:105) L (Ω) + (cid:104) y ∂ , u ∂ (cid:105) ∂ Ω . (20)In practice, this becomes explicit with the constitutive relation G ( f , e ) = 0 as it will be seen inSection 3.4 (and more generally in [46, 47]). Note that this latter relation has to be accurately discretizedto ensure that the discretized power balance mimics the continuous one. Remark 5
By adding resistive port(s), dissipation(s) can easily be taken into account (both internal orat the boundary), as done in the finite-dimensional case via R ( x ) playing the role of a output feedbackgain matrix. In this case, the system becomes a parabolic system , the dissipative constitutive relation being represented by G . See [48, 49] for a detailed discussion about structure-preserving discretization ofdissipative systems. The Partitioned Finite Element Method (PFEM)
We are now in a position to introduce a general methodology to discretize infinite-dimensional pHs ina structure-preserving manner. The main contribution in this section is the application of PFEM to ageneral abstract class of pHs, unifying the previously published results. This generality is notably ofparticular interest for the development of a well-structured software for the numerical simulations ofphysics-based models. The power balances (14), (20) are deeply linked to a linear underlying Stokes-Dirac. The main idea of PFEM is to mimic this structure, in order to obtain a discretized copy ofthese power balances as (2). This systematically translates the Stokes-Dirac structure into a finite-dimensional Dirac structure. The compatible discretization, with respect to this Dirac structure, of theconstitutive relations allows to mimic the continuous power-balance. This method goes under the namePartitioned Finite Element Method (PFEM), and was originally presented in [19]. The procedure is anatural extension of MFEM to pHs and boils down to these three simple steps:1. System (5) is written in weak form;2. The integration by parts (7) is carried out on a partition of the system (5) to make the appropriateboundary control appear;3. A Galerkin method is employed to obtain a finite-dimensional system. For the approximation basis,the finite element method is used here but spectral methods can be chosen as well.This strategy of structured discretization in order to mimic the continuous power balance at the dis-crete level has been addressed for closed abstract linear hyperbolic systems in [28]. This pioneering workalready proposed the key point of PFEM: the integration by parts on a partition of the weak formulationof the system. The author called the obtained systems primal–dual or dual–primal formulation, depend-ing on which line is integrated by parts. In the port-Hamiltonian formalism, systems are opened withcontrol and observation. It appears that [28] admits PFEM as a generalization for structure-preservingspace discretization. The choice of a control in the pHs community is called a causality , and primal–dual or dual–primal correspond in this work to the canonical causalities (10) and (11) respectively. Consider smooth test functions v and v and the weak form of (5): (cid:104) v , f (cid:105) L (Ω , A ) = (cid:104) v , −L ∗ e (cid:105) L (Ω , A ) , (cid:104) v , f (cid:105) L (Ω , B ) = (cid:104) v , L e (cid:105) L (Ω , B ) . (21)Next the integration by parts is performed either on the first or on the second line (the system is partitioned ), depending on the causality. Integration by parts of the term (cid:104) v , −L ∗ e (cid:105) L (Ω , A ) In this case case, using (7), it is obtained: − (cid:104) v , L ∗ e (cid:105) L (Ω , A ) = − (cid:104)L v , e (cid:105) L (Ω , B ) + (cid:104) Γ v , Γ ⊥ e (cid:105) ∂ Ω . The boundary variable e ∂ := Γ ⊥ e in (10) explicitly appears. Then the equation defining the corre-sponding f ∂ := Γ e is put into weak form to obtain the final system for all smooth test functions v , v , and v ∂ : (cid:104) v , f (cid:105) L (Ω , A ) = − (cid:104)L v , e (cid:105) L (Ω , B ) + (cid:104) Γ v , e ∂ (cid:105) ∂ Ω , (cid:104) v , f (cid:105) L (Ω , B ) = (cid:104) v , L e (cid:105) L (Ω , B ) , (cid:104) v ∂ , f ∂ (cid:105) ∂ Ω = − (cid:104) v ∂ , Γ e (cid:105) ∂ Ω . (22)8ow, a Galerkin discretization is introduced. Test, energy and co-energy functions with the same sub-script are discretized using the same basis, for all t ≥ : (cid:3)(cid:3)(cid:3) ( t, x ) ≈ (cid:3)(cid:3)(cid:3) d ( t, x ) := N (cid:88) i =1 ϕ i ( x ) (cid:3) i ( t ) , ∀ x ∈ Ω , (cid:3)(cid:3)(cid:3) ( t, x ) ≈ (cid:3)(cid:3)(cid:3) d ( t, x ) := N (cid:88) i =1 ϕ i ( x ) (cid:3) i ( t ) , ∀ x ∈ Ω , (cid:3)(cid:3)(cid:3) ∂ ( t, s ) ≈ (cid:3)(cid:3)(cid:3) d∂ ( t, s ) := N ∂ (cid:88) i =1 ψ i∂ ( s ) (cid:3) i∂ ( t ) , ∀ s ∈ ∂ Ω , (23)where (cid:3) stands for v , f , and e and ϕ i ∈ H L , ϕ i ∈ L (Ω , B ) , and ψ i∂ ∈ L ( ∂ Ω , R m ) . Remark 6
In general, a discretization in the same basis of either ( f , e ) and ( f , e ) or ( f , f ) and ( e , e ) (as done in [33]) must be performed. The former is our choice since it directly leads to square massmatrices, while the latter may be more appropriate when dealing for instance with Maxwell’s equationsfor electromagnetics, see [42] and references therein for details on the difficulties that may then occur. Then plugging the approximations into (22), it is computed: M
00 0 M ∂ f f f ∂ = − D (cid:62)L B ⊥ D L − B (cid:62)⊥ e e e ∂ , (24)where vectors f , f , e , e , f ∂ , and e ∂ are given by the column-wise concatenation of the respectivedegrees of freedom of f d , f d , e d , e d , f d∂ , and e d∂ , and where the matrices are defined as follows: M ij = (cid:68) ϕ i , ϕ j (cid:69) L (Ω , A ) ,M mn = (cid:104) ϕ m , ϕ n (cid:105) L (Ω , B ) ,M lk∂ = (cid:10) ψ l∂ , ψ k∂ (cid:11) ∂ Ω , D mi L = (cid:10) ϕ m , L ϕ i (cid:11) L (Ω , B ) ,B ik ⊥ = (cid:10) Γ ϕ i , ψ k∂ (cid:11) ∂ Ω , (25)where ≤ i, j ≤ N , ≤ m, n ≤ N , and ≤ l, k ≤ N ∂ . System (24) is a kernel representation of aDirac structure as in (3) (see Remark 1). Remark 7
Note that matrices D L and B ⊥ are not square. The discrete Hamiltonian is naturally defined as the continuous one evaluated in the discrete energyvariables. As done in Section 2, it is easier to distinguish the linear hyperbolic from the parabolic case.
Hyperbolic case
In this setting, the flows f i , i = 1 , , are given by the time derivative of theenergy variables α i . Hence, the discretization of these energy variables is given by: α di ( t, x ) = α di (0 , x ) + (cid:90) t f di ( s, x ) d s, i = 1 , . The discrete Hamiltonian is then defined by H d ( α , α ) := H ( α d , α d ) , where α and α are the column-wise concatenation of the time varying coefficients of α d and α d in their respective basis. Definition 2
The discretization of the constitutive relations is said to be compatible if and only if: ∇ H d = (cid:18) ∇ α H d ∇ α H d (cid:19) = (cid:20) M
00 M (cid:21) (cid:18) e e (cid:19) . Proposition 1
If the discretisation of the constitutive relations is compatible, the discrete power balancereads at the discrete level: ddt H d = − ( e ∂ ) (cid:62) M ∂ f ∂ , which perfectly mimics the continuous identity. roof 1 A straightforward computation gives: ddt H d = ( ∇ α H d ) (cid:62) ddt α + ( ∇ α H d ) (cid:62) ddt α , = ( M e ) (cid:62) f + ( M e ) (cid:62) f , = ( e ) (cid:62) M f + ( e ) (cid:62) M f , = − ( e ∂ ) (cid:62) M ∂ f ∂ , where the symmetry of the mass matrices and the Dirac structure have been used. Remark 8
In the special case of linear hyperbolic systems, it has been seen that the co-energy formula-tion allows to take the constitutive relations into account directly in the differential equations. ApplyingPFEM to (18) then leads to an ODE, and the constitutive relations are then automatically discretized ina compatible manner.
Parabolic case
In this setting, only the flow f is the time derivative of the energy variable α .This energy variable is discretized as in the hyperbolic case. The discrete Hamiltonian is then definedby H d ( α ) := H ( α d ) . Definition 3
The discretization of the constitutive relation is said to be compatible if and only if: ∇ H d = M e . Proposition 2
If the discretisation of the constitutive relations is compatible, the discrete power balancereads at the discrete level: ddt H d = − e (cid:62) M f − e (cid:62) ∂ M ∂ f ∂ , which perfectly mimics the continuous identity. Proof 2
The proof can be derived similarly as in the hyperbolic case.
Remark 9
Of course, an accurate discretization of the implicit constitutive relation G ( f , e ) = 0 isalso required to conclude. This will be illustrated in Section 3.4. Integration by parts of the term (cid:104) v , L e (cid:105) L (Ω , B ) Using (7), it comes: (cid:104) v , L e (cid:105) L (Ω , B ) = (cid:104)L ∗ v , e (cid:105) L (Ω , A ) + (cid:104) Γ ⊥ v , Γ e (cid:105) ∂ Ω . Now the boundary variable e ∂ := Γ e explicitly appears, i.e. the causality considered in (11). Theweak formulation then reads: (cid:104) v , f (cid:105) L (Ω , A ) = − (cid:104) v , L ∗ e (cid:105) L (Ω , A ) , (cid:104) v , f (cid:105) L (Ω , B ) = (cid:104)L ∗ v , e (cid:105) L (Ω , A ) − (cid:104) Γ ⊥ v , e ∂ (cid:105) ∂ Ω , (cid:104) v ∂ , f ∂ (cid:105) ∂ Ω = (cid:104) v ∂ , Γ ⊥ e (cid:105) ∂ Ω . (26)Plugging the approximations (23) into (26), this time with ϕ i ∈ L (Ω , A ) , ϕ i ∈ H −L ∗ , and ψ i∂ ∈ L ( ∂ Ω , R m ) , gives the following kernel representation of a finite-dimensional Dirac structure: M
00 0 M ∂ f f f ∂ = −L ∗ − D (cid:62)−L ∗ − B (cid:62) e e e ∂ , (27)where the matrices D −L ∗ and B are defined by: D im −L ∗ = (cid:10) ϕ i , −L ∗ ϕ m (cid:11) L (Ω , A ) , B mk = (cid:10) Γ ⊥ ϕ m , ψ k∂ (cid:11) ∂ Ω . (28)The power balances proven above still hold true with this causality, where the role played by e ∂ and f ∂ have been switched.In the sequel, this methodology is applied to the wave equation, the Mindlin-Reissner plate model andthe heat equation. These models have been chosen to demonstrate the versatility of our methodology.The wave equation is the prototype of linear hyperbolic systems, and the first example treated byPFEM [19]. The Mindlin model combines wave dynamics and plane elastodynamics, and requires theintroduction of tensor-valued variables. Finally, the heat equation is the prototype of parabolic systems,and leads to a pHs with intrinsic algebraic constraint, namely, to a pHDAE.10 .2 The wave equation The wave equation is a well-known model, used as the first example of linear hyperbolic systems in manylecture notes and books. This work is no exception to the rule. However, to account for more realisticphysics, let us consider the heterogeneous and anisotropic multidimensional wave equation. The equationreads (see [35]): ρ ∂ w∂t = div( T grad( w )) + f, ( x , t ) ∈ Ω × [0 , t f ] , Ω ⊂ R N , (29)where ρ is the mass density (bounded from above and below), T is the tensorial Young modulus (sym-metric and positive definite) and w is the deflection from the equilibrium. The field f accounts fordistributed force, such as gravity.Let us denote α p := ρ ∂w∂t the linear momentum and α q := grad( w ) the strain, as energy variables.Hence the Hamiltonian is given as the total energy (summing kinetic and potential energies) by: H = 12 (cid:90) Ω (cid:40) α p ρ + α q · ( T α q ) (cid:41) dΩ . (30)The co-energy variables are by definition the variational derivatives of H with respect to the energyvariables, i.e. : e p := δ α p H = α p ρ = ∂w∂t , e q := δ α q H = T α q = T grad( w ) , (31)the velocity and stress. With these notations, equation (29) rewrites: ∂∂t (cid:18) α p α q (cid:19) = (cid:20) (cid:21) (cid:18) e p e q (cid:19) , (32)together with the constitutive relations given in (31).Let us denote: e := e p , e := e q , M := ρ, M := T − , L := grad , Γ ⊥ := γ ⊥ = γ · n , Γ := γ , where γ is the Dirichlet trace operator. Then the Hamiltonian (30) rewrites as (17). The wave equa-tion (29) with Neumann boundary control u ∂ = γ ⊥ ( e q ) is given by (18).The application of PFEM directly gives: (cid:20) M ρ M T − (cid:21) ddt (cid:18) e e (cid:19) = (cid:20) − D (cid:62) grad D grad (cid:21) (cid:18) e e (cid:19) + (cid:20) B ⊥ (cid:21) u ∂ M ∂ y ∂ = (cid:2) B (cid:62)⊥ (cid:3) (cid:18) e e (cid:19) , (33)where: M ijρ := (cid:68) ϕ i , ρ ϕ j (cid:69) L (Ω , R ) , M kl T − := (cid:10) ϕ k , T − ϕ l (cid:11) L (Ω , R N ) , are the discretizations of the operators M and M respectively. The Mindlin model is a generalization to the 2D case of the Timoshenko beam model and is expressedby a system of two coupled PDEs (see [50]): ρb ∂ w∂t = div( q ) + f, ( x , t ) ∈ Ω × [0 , t f ] , Ω ⊂ R ,ρb ∂ θ ∂t = q + Div( M ) + τ , (34)where ρ is the mass density, b the plate thickness, w the vertical displacement, θ = ( θ x , θ y ) (cid:62) collects thedeflection of the cross section along axes x and y respectively. The fields f, τ represent distributed forces11nd torques. Variables M , q represent the momenta tensor and the shear stress. Hooke’s law relatesthose to the curvature tensor and shear deformation vector: M := D b K ∈ S , q := D s γ , K := Grad( θ ) ∈ S , γ := grad( w ) − θ .D s := E Y bK sh ν ) is the shear rigidity coefficient, where E Y is the Young modulus, ν is the Poisson modulus, K sh is the shear correction factor. Tensor D b is the bending stiffness: D b ( · ) = E Y b − ν ) [(1 − ν )( · ) + ν Tr( · )] . (35)An appropriate selection of the energy variables is the following [14, 13]: α w = ρb∂ t w, α θ := ρb ∂ t θ , A κ := K , α γ := γ . (36)The Hamiltonian H (total energy) is expressed in terms of energy variables as: H = 12 (cid:90) Ω (cid:26) ρb α w + 12 ρb (cid:107) α θ (cid:107) + A κ .. ( D b A κ ) + D s (cid:107) α γ (cid:107) (cid:27) Ω , (37)where A .. B denotes the tensor contraction. The co-energy variables are defined as follows: e w := δ α w H = ∂ t w, e θ := δ α θ H = ∂ t θ , E κ := δ A κ H = M , e γ := δ α γ H = q . (38)System (34) is then expressed in port-Hamiltonian form as [14] (forces and torques have been omittedfor simplicity): ∂∂t α w α θ A κ α γ = Div I × Grad grad − I × e w e θ E κ e γ . (39)By applying the divergence theorem, the energy rate is expressed as the duality product of the boundaryvariables: dHdt = (cid:104) ∂ t α w , e w (cid:105) L (Ω , R ) + (cid:104) ∂ t α θ , e θ (cid:105) L (Ω , R ) + (cid:104) ∂ t A κ , E κ (cid:105) L (Ω , S ) + (cid:104) ∂ t α γ , e γ (cid:105) L (Ω , R ) , = (cid:104) γ e w , γ ⊥ e γ (cid:105) ∂ Ω + (cid:104) γ e θ , γ ⊥ E κ (cid:105) ∂ Ω = (cid:104) y ∂, , u ∂, (cid:105) ∂ Ω + (cid:104) y ∂, , u ∂, (cid:105) ∂ Ω , (40)where: u ∂, = γ ⊥ e γ , u ∂, = γ ⊥ E κ , y ∂, = γ e w , y ∂, = γ e θ . The traces γ u = u | ∂ Ω , γ ⊥ U = U · n | ∂ Ω correspond to the Dirichlet trace for R d vectors and to thenormal trace for R d × d tensors. The mass operators are given by: M = (cid:20) ρb ρb (cid:21) , M = (cid:20) D − b D − s (cid:21) . (41)The L , Γ , and Γ ⊥ operators are: L = (cid:20) Gradgrad − I × (cid:21) , Γ = (cid:20) γ γ (cid:21) , Γ ⊥ = (cid:20) γ ⊥ γ ⊥ (cid:21) . (42)Introducing the approximations for the test and co-energy variables: (cid:52) w = N w (cid:88) i =1 φ iw (cid:52) iw , (cid:52) θ = N θ (cid:88) i =1 φ iθ (cid:52) iθ , (cid:52) κ = N κ (cid:88) i =1 Φ iκ (cid:52) iκ , (cid:52) γ = N γ (cid:88) i =1 φ iγ (cid:52) iγ , (43)where (cid:52) = { v, e } , and for the boundary controls: (cid:3) ∂, = N ∂, (cid:88) i =1 ψ i∂, (cid:3) i∂, , (cid:3)(cid:3)(cid:3) ∂, = N ∂, (cid:88) i =1 ψ i∂, (cid:3) i∂, , (cid:3) = { v, u, y } , ψ i∂, ∈ R , ψ i∂, ∈ R , (44)12FEM can be applied to obtain: Diag M ρh M I θ M D − b M D − s ˙ e w ˙ e θ ˙ e κ ˙ e γ = − D (cid:62) grad − D (cid:62) Grad − D (cid:62) Grad grad D e w e θ e κ e γ + B ⊥ ,w
00 B ⊥ ,θ (cid:18) u ∂, u ∂, (cid:19) , Diag (cid:20) M ∂, M ∂, (cid:21) (cid:18) y ∂, y ∂, (cid:19) = (cid:20) B (cid:62)⊥ ,w (cid:62)⊥ ,θ (cid:21) e w e θ e κ e γ . (45)The notation Diag denotes a block diagonal matrix. The mass matrices M ρh , M I θ , M C b , M C s arecomputed as: M ijρh = (cid:10) φ iw , ρhφ jw (cid:11) L (Ω) ,M mnI θ = (cid:104) φ mκ , I θ φ nκ (cid:105) L (Ω , R ) , M pq D − b = (cid:10) Φ pκ , D − b Φ qκ (cid:11) L (Ω , R × sym ) ,M rsD − s = (cid:10) φ rγ , D − s φ sγ (cid:11) L (Ω , R ) , (46)where i, j ∈ { , N w } , m, n ∈ { , N θ } , p, q ∈ { , N κ } , r, s ∈ { , N γ } . Matrices D grad , D Grad , D assumethe form: D rj grad = (cid:10) φ rγ , grad φ jw (cid:11) L (Ω , R ) ,D pn Grad = (cid:104) Φ pκ , Grad φ nθ (cid:105) L (Ω , R × sym ) , D rn = − (cid:10) φ rγ , φ nθ (cid:11) L (Ω , R ) . (47)Matrices B w , B θ , M ∂, , M ∂, are computed as: B ig ⊥ ,w = (cid:68) γ φ iw , ψ g∂, (cid:69) ∂ Ω , B mg ⊥ ,θ = (cid:68) γ φ mθ , ψ g∂, (cid:69) ∂ Ω , M fg∂, = (cid:68) ψ f∂, , ψ g∂, (cid:69) L ( ∂ Ω , R m ) ,M lk∂ = (cid:10) ψ l∂, , ψ k∂, (cid:11) L ( ∂ Ω , R m ) , f, g ∈ { , N ∂, } ,l, k ∈ { , N ∂, } (48)The discrete Hamiltonian is then computed as: H d = 12 e (cid:62) w M ρh e w + 12 e (cid:62) θ M I θ e θ + 12 e (cid:62) κ M D − b e κ + 12 e (cid:62) γ M D − s e γ . (49)From system (45) the discrete energy rate is readily obtained: dH d dt = y (cid:62) ∂, M ∂, u ∂, + y (cid:62) ∂, M ∂, u ∂, . (50)The discrete energy rate then mimics its infinite dimensional counterpart. Remark 10
Equivalently a purely mixed formulation can be obtained by integrating by parts the thirdand fourth lines of (39) . In this case, the system of equations gathers together a plane elasticity problem[6] and a wave equation in mixed form. Conforming finite elements for the plane elasticity system onsimplicial meshes have been constructed in [7]. The simpler PEERS elements based on a weak symmetryformulation have been proposed in [5]. The PEERS elements have been used in [10] to construct a stablelocking-free mixed formulation for the static Mindlin problem.
The heat equation is the simplest example of parabolic system. Instead of rewriting the well-known PDEunder a pHs, a direct pHs modelling is presented, as done in [46, 47]. The model is constructed in orderto keep apart thermodynamical principles from equations of state. Indeed, the pHs formalism allows tomodify the latter, by keeping the structure of the former.Let Ω ⊂ R N be a bounded open connected set. Assume that this domain models a rigid body:its volume does not change over time and no chemical reaction is to be found. Let us denote: ρ themass density, u the internal energy density, J Q the heat flux, T the local temperature, β := 1 T thereciprocal temperature, s the entropy density, J S := β J Q the entropy flux, C V := (cid:0) dudT (cid:1) V the isochoricheat capacity. 13he first law of thermodynamic reads: ρ ∂u∂t = − div( J Q ) . (51)Under the hypothesis of an inert rigid solid, Gibbs formula reads d u = T d s , giving: ∂u∂t = T ∂s∂t . (52)Defining σ := grad ( β ) J Q , and seeing u as a function of the entropy density s , Gibbs formula (52) gives: ρ ∂s∂t = − div( J S ) + σ. (53)Then σ is the irreversible entropy production.In this work, the following constitutive equations of state will be assumed:• The rigid body is at room temperature: the Dulong-Petit model is supposed to be satisfied, i.e. u = C V T , with time-invariant C V ;• The thermal conduction is given by Fourier’s law, with a symmetric positive tensor λ : J Q = − λ grad( T ) .Thanks to (51) and the equations of state, we easily recover the classical PDE for the temperature T : ρC V ∂T∂t = div( λ grad( T )) .The “ L -energy” (cid:0)(cid:82) Ω ρC V T dΩ (cid:1) is commonly used as Hamiltonian for the heat equation. However,it lacks of a thermodynamical meaning. The internal energy would be more accurate for this physicalproblem, even though it will rise some difficulties. Nevertheless, the pHs formalism allows dealing withit, and PFEM proves to be powerful enough to discretize the system in a structure-preserving mannereven for this choice of Hamiltonian.Let the internal energy be seen as a functional of the local entropy as energy variable: α s := ρs , then: H := (cid:90) Ω ρu ( α s )dΩ . The co-energy variable is given by e s := δ α s H = d ρu d ρs = T , the local temperature. Denoting e S := J S ,one can introduce a new flow variable f S such that: (cid:18) ∂α s ∂t f S (cid:19) = (cid:20) − div − grad 0 (cid:21) (cid:18) e s e S (cid:19) + (cid:18) σ (cid:19) . Obviously, f S = − grad( e s ) . In order to get a formally skew-symmetric operator, let us also introducean entropy port ( f σ , e σ ) , such that e σ = − σ . Then: ∂α s ∂t f S f σ = − div − − grad 0 01 0 0 e s e S e σ . (54) Remark 11
As surprising as it can be, in this setting, Fourier’s law appears to be stated in a nonlinear way: e s e S − λf S = 0 . This comes from the necessity to express the constitutive relations in function ofthe flows and efforts appearing in the equation defining the Stokes-Dirac structure. Remark 12
Two variables have been added to obtain (54) , but only one equation naturally appears: f σ = e s . Thus, another equation is needed to close the system: G ( f , e ) = 0 . Here, it is given bythe definition of the irreversible entropy production σ := grad( β ) J Q , rewritten in the flows and effortsvariables. This leads to the nonlinear constitutive relation: f S e S + f σ e σ = 0 . Remark 13
Usually the system energy is taken to be (cid:0)(cid:82) Ω ρC V T dΩ (cid:1) , that gives rise to the well-knownlinear diffusive system. At first glance, our approach may be surprising since it leads to a losslessnonlinear differential-algebraic system. There is indeed a major advantage in doing so, in view of themodelling and discretization of complex systems by interconnection of several pHs. For instance, if onewants to interconnect both thermal and mechanical processes, for the energy exchanges to be consistent,physics must be coherent. When dissipation occurs, through e.g. friction or viscosity, the kinetic energyis converted into internal energy. Hence, for a physically meaningful system, the internal energy is indeedthe one to consider. f := ∂α s ∂t , f := (cid:18) f S f σ (cid:19) , e := e s , e := (cid:18) e S e σ (cid:19) , and: L := (cid:18) − grad1 (cid:19) , Γ ⊥ := (cid:0) − γ ⊥ (cid:1) , Γ := γ . Then, the heat equation (54) with boundary control e ∂ := u ∂ = Γ e = γ ( e s ) and boundary observation f ∂ := − y ∂ = Γ ⊥ e = γ ⊥ ( e S ) rewrites under the form (11). Thus PFEM will be applied with an integration by parts on the second line in this strategy, leading to (27), which rewrites with the currentvariables: M s S σ
00 0 0 M ∂ ddt α s f S f σ − y ∂ = − div − M σ − D (cid:62)− div M σ − B (cid:62) e s e S e σ u ∂ , (55)where: M ijs := (cid:68) ϕ i , ϕ j (cid:69) L (Ω , R ) , M klS := (cid:10) ϕ kS , ϕ lS (cid:11) L (Ω , R N ) ,M (cid:101) k (cid:101) lσ := (cid:68) ϕ (cid:101) kσ , ϕ (cid:101) lσ (cid:69) L (Ω , R N ) , B km := − (cid:10) Γ ⊥ ϕ kS , ψ m∂ (cid:11) L ( ∂ Ω , R ) , with ϕ k := (cid:18) ϕ kS (cid:19) if ≤ k ≤ N S , and ϕ k := (cid:18) ϕ k − N S σ (cid:19) if N S + 1 ≤ k ≤ N S + N σ .To be compatible, the discretizations of the constitutive relations are given as follows:• Dulong-Petit model reads: M s α s = M ρC V e s , with M ijρC V := (cid:68) ϕ i , ρ C V ϕ j (cid:69) L (Ω , R ) ;• Following Remark 11, Fourier’s law reads: Λ f S = M e s e S , with Λ ij := (cid:68) ϕ iS , λ ϕ jS (cid:69) L (Ω , R N ) and M ije s := (cid:68) ϕ iS , e s ϕ jS (cid:69) L (Ω , R N ) ;• The constitutive law coming from the introduction of the irreversible entropy production, as ex-plained in Remark 12, is taken into account by: ( e S ) (cid:62) M S f S + ( e σ ) (cid:62) M σ f σ = 0 . (56) Remark 14
In Fourier’s law, the mass matrix M e s depends on the co-energy variable e s . This will risedifficulties for the numerical solution in time. To conclude, the structure-preserving property can be appreciated in the following result.
Proposition 3
Let H d ( α s ) := H ( α ds ) be the discrete Hamiltonian, where α ds is the discretization of theenergy variable in the basis ϕ . It holds: ddt H d = u (cid:62) ∂ M ∂ y ∂ , that is the first law of thermodynamics at the discrete level. Proof 3
Thanks to the compatible discretization of the Dulong-Petit model, Proposition 2 gives: ddt H d = − e (cid:62) M f − e (cid:62) ∂ M ∂ f ∂ . By definition of f , e , M , e ∂ , and f ∂ one computes: ddt H d = − e (cid:62) M f − e (cid:62) ∂ M ∂ f ∂ , = − (cid:18) e S e σ (cid:19) (cid:62) (cid:20) M S
00 M σ (cid:21) (cid:18) f S f σ (cid:19) + u (cid:62) ∂ M ∂ y ∂ , = u (cid:62) ∂ M ∂ y ∂ , thanks to the constitutive relation (56) coming from the irreversible entropy production. emark 15 Fourier’s law does not contribute to the power balance of the internal energy. Nevertheless,such a constitutive relation is needed for the problem to be well-defined.
Remark 16
The methodology detailed so far is certainly not limited to the previous three examples.Indeed higher-order differential [15], curl operator for Maxwell’s equations [42], nonlinear system [20],and different Hamiltonian choices can be handled as well. For instance, in the case of the heat equation,the entropy or the classical L -norm of the temperature can be alternatively considered as Hamiltonianfunctional [46, 47]. In addition, mixed boundary conditions can be incorporated either by introducingLagrange multipliers or by employing a virtual domain decomposition method [16]. As already mentioned,dissipation (both in the domain and on the boundary) can also be considered in this strategy [48, 49].Hence a very wide class of (nonlinear) multiphysics systems can be discretized in a structure-preservingmanner (with well-represented exchanges of energy between the subsystems). In the next section we present an ongoing project which has been initiated to prove the efficiencyof the PFEM methodology, leveraging well-established and robust software tools for the finite elementdiscretization of partial differential equations and time integration.
In this section the main features related to the numerical simulation of pHs in the framework of theongoing project named
SCRIMP (Simulation and ContRol of Interactions in Multi-Physics) are detailed.The aim is to provide a flexible prototype
Python code for the numerical simulation of pHs both forresearch and educational purposes. In addition to numerical experiments proposed later in Section 5,the reader is referred to interactive companion
Jupyter notebooks [17] to learn how to numericallysolve the model problems introduced in Section 3 with
SCRIMP . In the following, the key ideas behind
SCRIMP are mentioned and then a specific emphasis on both space and time discretizations is given.
SCRIMP
In short, the key ideas related to the design of
SCRIMP are provided:• The
Python dynamic programming language has been selected due to its expressiveness and theavailability of high-level interfaces to scientific computing software libraries [37];•
SCRIMP assumes to rely on open-source, external software for the finite dimensional discretizationof partial differential equations;•
SCRIMP encapsulates the finite dimensional objects related to the finite element discretization inspace (e.g. matrices) to deduce the resulting linear or nonlinear pHs in a generic pHODE/pHDAEform as proposed in [9];• For multiphysics problems, this design offers the advantage that discretization in space may behandled by different software components depending on the discipline or on the modelling. Themodularity and the object-oriented nature of
Python thus offer the flexibility to easily combinethe different pHs to deduce the global interconnected system. This is much in line with themathematical theory of pHs [54]. Furthermore we note that interconnections of different systems(with e.g. the transformer or gyrator transformations [54]) can be easily incorporated.The design of
SCRIMP is based on procedural and object-oriented paradigms and thus follows thestandard ideas governing most of the numerical PDE software. Whereas a detailed exposition of thedesign patterns of
SCRIMP and its performance will be published elsewhere, concrete illustrations ofmost of these key ideas can be found in the companion
Jupyter notebooks [17]. The description of thecurrent numerical methods related to space and time discretizations available in
SCRIMP is given.
As outlined in Section 3, PFEM relies on an abstract variational formulation written in appropriate finiteelement spaces. 16o perform the semi-discretization in space, we rely on
FEniCS [2], an open-source
C++ scientificsoftware library that provides a high-level
Python interface. The
FEniCS
Project is mainly basedon a collection of software components targeting the automated solution of partial differential equa-tions via the finite element method. Its core components notably include the Unified Form Language(
UFL ) [3], the
FEniCS
Form Compiler (
FFC ) [30] and the finite element library
DOLFIN [39], whichcontains various types of conforming finite element methods, e.g., nodal Lagrangian finite elements forgrad-conforming approximations or non non-nodal finite elements (e.g., Raviart-Thomas spaces for div-conforming approximations) as well. These families of finite elements are notably required to tackle thediscretization in space of our core problems.A key point to facilitate the generic implementation of PFEM is the use of
UFL . UFL is indeed anexpressive domain-specific language for abstractly representing (finite element) variational formulationsof differential equations. In particular, this language defines a syntax for the integration of variationalforms over various domains. This simply leads to an expressive implementation that is close to theabstract mathematical formulations presented in Section 3. The
FEniCS
Form Compiler
FFC thengenerates specialized
C++ code from the symbolic
UFL representation of variational forms and finiteelement spaces. The combination of these core elements makes
FEniCS a versatile and efficient softwarefor the finite element approximation of partial differential equations as outlined in [38]. Additionally,
FEniCS also provides an interface for state-of-the-art linear solvers and preconditioners from freelyavailable third-party libraries such as
PETSc [8]. This last feature may be especially useful to handlethe numerical simulation of large-scale pHs.
As outlined in Section 3, the semi-discretization in space of the resulting pHs leads to systems of ei-ther ordinary differential equations (ODE) or differential algebraic equations (DAE). Hence reliable andaccurate time integration methods must be provided.To offer a large panel of numerical methods, a high-level interface to well-established time integrationlibraries is provided in
SCRIMP . Concerning the numerical solution of ODEs, we provide light inter-faces to the
Assimulo library [4] and to the
SciPy time integration method scipy.integrate.solve_ivp that both include standard multistep and one-step methods for stiff and non-stiff ordinary differentialequations given in explicit form y (cid:48) = f ( t, y ) with y ( t ) = y where t and y denote the initial timeand initial condition, respectively. This formulation requires the solution of linear systems of equationsinvolving sparse finite element mass matrices. State-of-the-art sparse direct solvers based on Gaussianfactorization are used for that purpose. Through Assimulo , the popular
CVODE solver from Sun-dials [27] is also accessible. For nonstiff problems,
CVODE relies on the Adams-Moulton formulas,with the order varying between and . For stiff problems, CVODE includes schemes based on Back-ward Differentiation Formulas (BDFs) with order varying between and . In addition, we also proposestandalone implementations of symplectic time integration methods such as the second order accurateStormer-Verlet method.The interface to Assimulo also allows one to handle the numerical solution of linear DAEs throughthe use of the
Sundials IDA solver . IDA is a package for the solution of differential algebraic equationsystems written in the form F ( t, y, y (cid:48) ) = 0 with y ( t ) = y . The integration method in IDA is basedon variable-order, variable-coefficient BDF in fixed-leading-coefficient form, where the method ordervaries between and . We note that setting the initial conditions properly is of utmost importancefor a DAE solver. To do so, we rely on the IDA_YA_YDP_INIT method to find consistent initialconditions for the time integration. We refer the reader to [27, Section 2.3] for additional details on
IDA .As standalone methods, we have considered the second-order accurate Stormer-Verlet and fourth-orderaccurate Runge-Kutta (RK4) methods for the solution of linear DAEs.To the best of our knowledge, open-source libraries for the solution of general nonlinear differentialalgebraic equations with high-level
Python interfaces are not yet available. Hence a simple forwardin time integration method for the solution of the nonlinear pHDAE related to the energy formulationof the heat equation problem has been provided; see [45] for illustrations and discussion. As a futuredirection, we plan to investigate the potential of the
PETSc ’s time stepping library TS [1] to be ableto tackle the solution of large-scale pHDAE systems. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html https://computing.llnl.gov/sites/default/files/public/cv_guide.pdf https://computing.llnl.gov/projects/sundials/ida .4 Model reduction of port-Hamiltonian systems Structure-preserving model reduction is of significant importance for stability analysis, optimization orcontrol of problems related to pHs. Hence structure-preserving model reduction algorithms have beenimplemented in
SCRIMP . In particular, the structure-preserving model reduction algorithm (Algorithm1) proposed in [21] has been selected in the pHODE case. We refer the reader to [17] for an illustration,where the model reduction of the pHs related to the wave equation problem is considered. While for linearpHDAE systems consolidated methodologies have been proposed (see, e.g., [24]), structure-preservingmodel reduction for general nonlinear differential algebraic systems remains to be explored, to the bestof our knowledge. This is a significant research direction to be considered within
SCRIMP in a nearfuture.
In this section, PFEM is applied to the pHs presented in Section 3. We specifically learn how to defineand solve those problems with
SCRIMP . These tutorials introduce the methodology step-by-step andare supposed to be self-contained and independent from the others. We refer the reader to the companion
Jupyter notebooks [17] for additional information.
We first recall the continuous problem related to the anisotropic heterogeneous wave equation, enrichedwith internal and boundary damping, and tackle the semi-discretization in space of the port-Hamiltoniansystem through the PFEM methodology. This discretization leads to a pHODE formulation as explainedin Section 3.2. After time discretization, we perform a numerical simulation to obtain an approximationof the space-time solution.
We consider the two-dimensional heterogeneous anisotropic wave equation with impedance boundarycondition defined for all ( t ≥ as: ρ ( x ) ∂ ∂t w ( t, x ) − div (cid:16) T ( x ) · grad w ( t, x ) (cid:17) = − (cid:15) ( x ) ∂ t w ( t, x ) , x ∈ Ω ,Z ( x ) ( T ( x ) · grad w ( t, x )) · n + ∂ t w ( t, x ) = 0 , x ∈ ∂ Ω ,w (0 , x ) = w ( x ) , x ∈ Ω , t = 0 ∂ t w (0 , x ) = w ( x ) , x ∈ Ω , t = 0 , with Ω ⊂ R an open bounded spatial domain with Lipschitz-continuous boundary ∂ Ω . We considerhere a rectangular shaped domain for Ω . w ( t, x ) denotes the deflection from the equilibrium positionat point x ∈ Ω and time t . ρ ∈ L ∞ (Ω) (positive and bounded from below) denotes the mass density, T ∈ L ∞ (Ω) × (symmetric and coercive) the Young’s elasticity modulus, (cid:15) a positive viscous dampingparameter and Z ( x ) is the positive impedance function defined on ∂ Ω , respectively. We initialize here the Python object related to the
Wave_2D class of
SCRIMP . This object will be usedthroughout this section.
W = SCRIMP.Wave_2D()
We define the constants related to the rectangular domain Ω . The coordinates of the bottom left ( x , y )and top right ( x L , y L ) corners of the rectangle are required. x0, xL, y0, yL = 0., 2., 0., 1. We then define the time interval related to the time discretization. t i , t f denote the initial and final timeinstants respectively. 18 i, tf = 0., 5. We specify that we choose the
Assimulo external library to be used later for the time integration of theresulting ODE and provide the value of the time step. This should be considered as a reference valuesince adaptative methods in time can be used later. dt = 1.e-3ode_library = 'ODE:Assimulo'
For the finite element discretization of the pHs, the
FEniCS library is used in the
Wave_2D class of
SCRIMP . Hence to properly use
FEniCS expression definition, we provide the definition of the differentvariables in
C++ code given in strings. We first specify the mass density as a function depending onthe space coordinates. Hence in this expression, x [0] corresponds to the first spatial variable and x [1] tothe second one, respectively. Rho = 'x[0]*x[0] * (2.-x[0])+ 1'
We specify the Young’s elasticity modulus tensor. Three components are only required due to thesymmetry property of this tensor.
T11 = 'x[0]*x[0]+1'T12 = 'x[1]'T22 = 'x[0]+2'
We finally set the impedance function Z defined on the boundary of the domain. Here a constant valueis used on ∂ Ω . We also provide the viscous damping parameter ( eps ). Z = '0.1'eps = '25 * x[0] * (xL-x[0]) * x[1] * (yL-x[1])'
Finally we specify the initial conditions of the problem related to the energy variables and to the deflec-tion.
Aq_0_1= '0'Aq_0_2= '0'Ap_0 = '0'W_0 = '0'
We are now able to completely define the problem at the continuous level. We start by specifying thatthe computational domain Ω is of rectangular shape. To do so, we provide the coordinates of the bottomleft and top right corners to the Wave_2D object.
W.Set_Rectangular_Domain(x0, xL, y0, yL);
Remark 17
General
Gmsh meshes can be imported by the user. However, for the time being, the librarydoes not allow the treatment of mixed boundary conditions on generic meshes.
We provide next the time integration interval.
W.Set_Initial_Final_Time(ti, tf);
We then provide the physical parameters related to the wave equation: the mass density, the Young’selasticity modulus tensor and the impendance function, respectively.19 .Set_Physical_Parameters(Rho, T11, T12, T22);
We then specify the complete modelling for the damping and thus provide information related to theimpedance function and viscous damping parameter, respectively.
W.Set_Damping(damp=['impedance', 'fluid'], Z=Z, eps=eps);
The user has to provide the temporal and spatial parts of the boundary control function (
Ub_tm0 and
Ub_sp0 , respectively).
W.Set_Boundary_Control(Ub_tm0=lambda t: np.sin( 2 * 2*pi/tf *t) * 25 ,␣ (cid:44) → Ub_sp0='x[0] * x[1] * (1-x[1])');
Finally we provide the initial conditions for the ODE.
W.Set_Initial_Data(Aq_0_1='0', Aq_0_2='0', Ap_0='0', W_0='0');
We start by selecting the computational mesh which is generated with
Gmsh and saved as a .xml file.Here the parameter rf n _ num corresponds to a mesh refinement parameter. W.Set_Gmsh_Mesh('rectangle.xml', rfn_num=2);
To perform the discretization in space, we must first specify the conforming finite element approximationspaces to be used (see [26]). Concerning the energy variables associated with the strain, we select theRaviart-Thomas finite element family known as RT k consisting of vector functions with a continuousnormal component across the interfaces between the elements of a mesh. For the energy variablesassociated with the linear momentum and the boundary variables, we choose the classical P k finiteelement approximation. The given combination of parameters rt_order=0, p_order=1, b_order=1 corresponds to the RT P P family. W.Set_Finite_Element_Spaces(family_q='RT', family_p='P', family_b='P',rq=0, rp=1,␣ (cid:44) → rb=1); We then perform the semi-discretization in space of the weak formulation with PFEM. At the endof this stage, the complete formulation of the pHODE is obtained. The different matrices related tothe pHODE system are constructed in the Assembly method of the
Wave_2D class of
SCRIMP andare directly accessible through the object of the
Wave_2D class. The finite element assembly relies onthe variational formulation of PFEM and exploits the level of abstraction provided by the UFL usedin
FEniCS , leading to a code that is close to the mathematical formulation. The divergence basedformulation is selected leading to a pHODE system. In other words, the integration by parts will beperformed on the second line of (32).
W.Assembly(formulation='Div');
To perform the time integration of the pHODE, we first need to interpolate both the control function onthe boundary and the initial data on the appropriate finite element spaces.
W.Project_Boundary_Control()W.Project_Initial_Data();
Then we specify the parameters related to the time discretization.
W.Set_Time_Setting(dt); https://gmsh.info/ .1.7 Numerical approximation of the space-time solution We are now able to perform the time integration of the resulting pHODE system and deduce the be-haviour of both the energy variables and the Hamiltonian with respect to the time and space variables,respectively. Detailed information from the
Assimulo library is included after time integration.
A, Hamiltonian = W.Time_Integration(ode_library)ODE Integration using assimulo built-in functions:Final Run Statistics: ---Number of steps : 614Number of function evaluations : 800Number of Jacobian*vector evaluations : 2977Number of function eval. due to Jacobian eval. : 0Number of error test failures : 0Number of nonlinear iterations : 797Number of nonlinear convergence failures : 51Solver options:Solver : CVodeLinear multistep method : BDFNonlinear solver : NewtonLinear solver type : SPGMRMaximal order : 3Tolerances (absolute) : 1e-05Tolerances (relative) : 1e-05Simulation interval : 0.0 - 5.0 seconds.Elapsed simulation time: 0.9727537930002654 seconds.
We represent the two-dimensional mesh with corresponding degrees of freedom for each variable in Figure1.
W.notebook = TrueW.Plot_Mesh_with_DOFs()
Figure 1: Two-dimensional triangular mesh with corresponding degrees of freedom for each variable forthe anisotropic wave problem with damping.We plot the Hamiltonian function versus time in Figure 2. Here tspan represents the collection ofdiscrete times due to the possibly adaptative time procedure used in the time integration library.21 .Plot_Hamiltonian(W.tspan,Hamiltonian, marker='o')
Figure 2: Hamiltonian function versus time for the anisotropic, heterogeneous and boundary controlledwave problem with damping.The behaviour of the deflection is graphically represented at a given time. Here we simply plot thedeflection at the final time of the simulation in Figure 3.
W.Plot_3D(W.Get_Deflection(A), tf, 'Deflection at t='+str(tf))
Figure 3: Deflection at the final time for the anisotropic, heterogeneous and boundary controlled waveproblem with damping.The related
Jupyter notebook [17] further illustrates how to obtain a structure-preserving reducedmodel of this port-Hamiltonian system. After application of the model reduction algorithm proposedin [21], a pHODE of reduced size has to be integrated to obtain an approximate solution of the wavepropagation problem. This is further illustrated on the simple application detailed in this section. Inaddition, a supplementary notebook illustrates the numerical simulation of the wave equation problem,when mixed boundary conditions (i.e. Dirichlet and Neumann conditions) on the boundary controlfunction are imposed by Lagrange multipliers [16].22 .2 The Mindlin plate problem
We first recall the considered continuous problem related to the Mindlin plate and tackle the semi-discretization in space of the pHs by PFEM. After transformation and time discretization, we performa numerical simulation to obtain an approximation of the space-time solution. As in Section 5.1, theprocedure is described step-by-step and detailed explanations and numerical illustrations are provided.
Consider the Mindlin plate problem defined for all t ≥ as: ρb∂ tt w = div( q ) , x ∈ Ω = { , L x } × { , L y } ,ρb ∂ tt θ = q + Div( M ) (57)with initial conditions: w (0 , x ) = 0 ,∂ t w (0 , x ) = w ( x ) , θ (0 , x ) = ,∂ t θ (0 , x ) = , x ∈ Ω , t = 0 , x ∈ Ω , t = 0 , (58)and boundary conditions: w | Γ D = u D ( t ) , q · n | Γ N = u N ( t ) , θ | Γ D = 0 , M · n | Γ N = 0 , Γ D = { x = 0 } , Γ N = { x = L x , y = 0 , y = L y } . (59)Mixed boundary conditions are considered in this example. The subsets Γ N , Γ D represent the subsets ofthe boundary where Neumann and Dirichlet conditions hold respectively. The Dirichlet conditions areenforced using Lagrange multipliers. The PFEM discretization then leads to a pHDAE, as explained in[16] for the wave equation. The following expressions have been considered for the initial and boundaryconditions: w = xy, u D = 0 . (cid:18) π tt f (cid:19) − , u N = 10 sin (cid:18) π xL x (cid:19) (cid:16) − exp − ttf (cid:17) . (60) We initialize here the
Python object related to the
Mindlin class of
SCRIMP . This object will be usedthroughout this section.
Min = SCRIMP.Mindlin()
We define the constants related to the rectangular domain Ω . The coordinates of the bottom left ( x , y ) = (0 , and the top right ( x L , y L ) = ( L x , L y ) corners of the rectangle are required. x0, xL, y0, yL = 0., 2., 0., 1. As in the previous example, the time interval related to the time discretization is defined as follows: ti, tf = 0., 0.01
A Runge-Kutta method for the time integration of the system is prescribed. This method is conditionallystable, so the time-step has to be set accurately to avoid numerical instabilities. dt = 1.e-6dae_library = 'DAE:RK4_Augmented'
FEniCS expressions definition
The
FEniCS library is also used in the
Mindlin class of
SCRIMP . The coefficients related to the physicalparameters of the isotropic plate can be provided as either real numbers or
FEniCS expressions.23 = 7*10**10rho = 2700nu = 0.3h = 0.1k = 5/6
Similarly the initial vertical condition w can be defined as a string. It represents a C++ code that willbe compiled by the
Dolfin library of
FEniCS . ew_0 = 'x[0]*x[1]' This means that the initial velocity satisfies w = xy . Note that the initial condition has to be compatiblewith the boundary conditions. The other initial conditions will be set to zero. We are now able to completely define the problem at the continuous level. We start by specifying thatthe computational domain Ω is of rectangular shape. To define Ω , we provide the coordinates of thebottom left and top right corners to the Mindlin object.
Min.Set_Rectangular_Domain(x0, xL, y0, yL);
The time integration interval is then given.
Min.Set_Initial_Final_Time(ti, tf);
The physical parameters related to the Mindlin plate are set.
Min.Set_Physical_Parameters(rho, h, E, nu, k, init_by_value=True);
Finally the initial conditions in terms of co-energy variables are also set.
Min.Set_Initial_Data(W_0='0', Th1_0='0', Th2_0='0',\ew_0=ew_0, eth1_0='0', eth2_0='0',\Ekap11_0='0', Ekap12_0='0', Ekap22_0='0',\egam1_0='0', egam2_0='0');
We start by selecting the computational mesh which is generated with
FEniCS inner mesh utilities. Thefirst parameter corresponds to a mesh refinement parameter.
Min.Generate_Mesh(10, structured_mesh=True);
To perform the discretization in space, the conforming finite element approximation spaces to be usedhas to be specified. The finite element for the linear and angular velocity are Lagrange polynomials oforder r . The momenta tensor and shear stress are chosen as Discontinuous Galerkin elements of order r − [13]. This choice of finite elements is similar to the one proposed in [22], but a simplicial mesh isused instead of a quadrilateral one. By default, the boundary variables are approximated as Lagrangepolynomials of order . This can be easily changed through the option family_b for the family, and rb for the degree. Min.Set_Finite_Elements_Spaces(r=1);
We then perform the semi-discretization in space of the weak formulation with PFEM. At the end ofthis stage, the complete formulation of the pHDAE is obtained. The different matrices related to thepHDAE system are constructed in the
Assembly_Mixed_BC method of the Mindlin class of
SCRIMP and are directly accessible through the object of the Mindlin class. The subsets named G1 , G2 , G3 , G4 ,denote the left, bottom, right and top sides of the rectangle, respectively.In SCRIMP the boundary control u ∂ is assumed to take the form:24 b_tm0(t) * Ub_sp0(x) + Ub_tm1(t) + Ub_sp1(x) Its derivative ˙ u ∂ is expressed as: Ub_tm0_dir(t) * Ub_sp0(x) + Ub_tm1_dir(t)
To integrate in time we need to provide the derivative of the boundary condition. This information isprovided by the variables
Ub_tm0_dir , Ub_tm1_dir , respectively. This is needed to reduce the resultingDAE of index to index . Min.Set_Mixed_Boundaries(Dir=['G1'], Nor=['G2', 'G3', 'G4'])Min.Assembly_Mixed_BC()Min.Set_Mixed_BC_Normal(Ub_tm0=lambda t: np.array([(1 - np.exp(-t/tf) ),0,0]) ,\Ub_sp0=('100000*sin(2*pi/xL*x[0])', '0.', '0.'))amp = 0.01omega = 2*pi/tfMin.Set_Mixed_BC_Dirichlet(Ub_tm0=lambda t : np.array([-amp*omega*np. (cid:44) → sin(omega*t),0,0]), Ub_sp0=('1.', '0.','0.'),\Ub_tm0_dir=lambda t : np.array([-amp*omega**2*np. (cid:44) → cos(omega*t),0,0])); To perform the time integration of the pHDAE, we first need to interpolate the boundary control functionand the initial data on the appropriate finite element spaces.
Min.Project_Boundary_Control()Min.Project_Initial_Data();
Finally the specification of the parameters related to the time discretization is made.
Min.Set_Time_Setting(dt);
For the numerical approximation of the solution of the pHDAE system, the algebraic condition isdifferentiated. The integrator 'DAE:RK4_Augmented' takes as input a pHDAE. Then, it exploits aprojection method to express the Lagrange multiplier in terms of the unknown [11], thus reducing theoriginal DAE system into a purely ODE one. This allows employing standard ODE solvers for the timeintegration, as discussed in Section 5.2.3.
A, Hamiltonian = Min.Time_Integration(dae_library)
Post-processing is performed similarly as in Section 5.1.8. Hence we omit the related
Python lines ofcode for sake of brevity. In Figure 4 the evolution of the Hamiltonian function is shown versus time.We note that the Dirichlet condition causes an increase in energy. In Figure 5 snapshots of the verticaldeflection at different instants are shown. We remark that the Neumann boundary condition causes theplate to bend asymmetrically.
This third tutorial aims at illustrating PFEM to discretize the pHs presented in Section 3.4, modellingthe heat equation. We specifically learn how to define and solve this problem with
SCRIMP . We firstdefine the continuous problem by using a specific class of
SCRIMP related to the heat equation in twodimensions. Then we tackle the discretization in space of the pHs through PFEM. The discretization ofthe energy formulation leads to a nonlinear pHDAE formulation. After time discretization, we performa numerical simulation to obtain an approximation of the space-time solution. Finally a simple post-processing is provided. 25igure 4: Hamiltonian versus time for the Mindlin plate problem. (a) Deflection w at t = t f / (b) Deflection w at t = t f / (c) Deflection w at t = 3 t f / (d) Deflection w at t = t f Figure 5: Snapshots of the vertical deflection of the Mindlin plate at different instants.26 .3.1 Problem statement
We consider the two-dimensional heterogeneous anisotropic heat equation defined for all t ≥ as ρ ( x ) C V ( x ) ∂∂t T ( t, x ) = div (cid:16) λ ( x ) · grad T ( t, x ) (cid:17) , x ∈ Ω ,T ( t, x ) = v ∂ ( t, x ) , x ∈ ∂ Ω ,T (0 , x ) = T ( x ) , x ∈ Ω , t = 0 with Ω ⊂ R an open bounded spatial domain with Lipschitz-continuous boundary ∂ Ω . v ∂ ( t, x ) representsthe boundary control function on the temperature (the notation u is kept for the internal energy density). T ( t, x ) denotes the temperature at point x ∈ Ω and time t . ρ ∈ L ∞ (Ω) (positive and bounded frombelow) denotes the mass density, λ ∈ L ∞ (Ω) × (symmetric and coercive) the thermal conductivity. Inthe following Ω is assumed to be of rectangular shape. We refer to [46, 47] for the modeling and discretization of various port-Hamiltonian formulations of thisproblem. The authors consider quadratic Lyapunov functional, entropy or internal energy as Hamilto-nian, respectively. We will consider the PFEM discretization of the internal energy functional formulationas proposed in Section 3.4, which will lead to a nonlinear pHDAE. Our goal in this tutorial is to showhow a pHDAE system can be formulated and solved with
SCRIMP . We initialize here the Python object related to the energy formulation of the
Heat_2D class of
SCRIMP ,that is assumed to be imported. This object will be used throughout this tutorial.
H = Energy()Energy corresponds to a class inherited from the
Heat_2D base class. This base class contains imple-mentations of the Lyapunov and entropy formulations as well.
The same lines of code as for the
Wave_2D and
Mindlin classes are used to define the constants relatedto the rectangular mesh. x0, xL, y0, yL = 0., 2., 0., 1.
The time interval related to the time discretization is specified similarly. ti, tf = 0., 5.
We provide the time step for the time discretization of the pHDAE as well. dt = 1.e-3
Using
FEniCS expressions, the physical parameters related to our model problem are defined. rho = 'x[0]*(xL-x[0]) + x[1]*(yL-x[1]) + 2'Lambda11 = '5. + x[0]*x[1]'Lambda12 = '(x[0]-x[1])*(x[0]-x[1])'Lambda22 = '3.+x[1]/(x[0]+1)'CV = '3.'
The initial conditions of the problem related to the temperature and to the flow and effort variables arethen given. The temperature follows a Gaussian behaviour for which we specify related parameters.27 mpl, sX, sY, X0, Y0 = 1000, xL/6, yL/6, xL/2, yL/2au_0 = '3000'eu_0 = ' ampl * exp(- pow( (x[0]-X0)/sX, 2) - pow( (x[1]-Y0)/sY, 2) ) + 1000'
The spatial part of the boundary control function is defined next.
Ub_sp0 = '''( abs(x[0]) <= DOLFIN_EPS ? 1. * (yL-x[1])*x[1] : 0 )+ ( abs(x[1]) <= DOLFIN_EPS ? 1. * (xL-x[0])*x[0] : 0 )+ ( abs(xL - x[0]) <= DOLFIN_EPS ? -15. * (yL-x[1])*x[1] : 0 )+ ( abs(yL - x[1]) <= DOLFIN_EPS ? 1. * (xL-x[0])*x[0] : 0 )'''
Finally we define the time-dependent part of the boundary control as a pure Python function. Thewhole boundary control function is then given as the product of the two quantities (
Ub_sp0 and
Ub_tm0 ,respectively). def Ub_tm0(t):if t<=2:return 500 * sin(2 * pi * t)else: return 0
We are now able to completely define the problem at the continuous level.
H.Set_Rectangular_Domain(x0, xL, y0, yL);H.Set_Initial_Final_Time(initial_time=ti, final_time=tf);H.Set_Physical_Parameters(rho=rho, Lambda11=Lambda11, Lambda12=Lambda12,␣ (cid:44) → Lambda22=Lambda22, CV=CV);
The structure-preserving discretization of the infinite-dimensional pHs with PFEM is described in detailin [47]. This leads to the pHDAE given in (55). The definition of the system at the discrete level followsthe same steps as for the two previous examples.
H.Set_Gmsh_Mesh(xmlfile='rectangle.xml', rfn_num=1);H.Set_Finite_Element_Spaces(family_q='RT', family_p='P', family_b='P',rq=0, rp=1,␣ (cid:44) → rb=1);H.Assembly(); To perform the time integration of the pHDAE, we first need to set and interpolate the initial data andthe boundary control function on the appropriate finite element spaces. Then, the time step is specified.
H.Set_Initial_Data(au_0=au_0, eu_0=eu_0, ampl=ampl, sX=sX, sY=sY, X0=X0, Y0=Y0);H.Project_Initial_Data();H.Set_Boundary_Control(Ub_tm0=Ub_tm0, Ub_sp0=Ub_sp0, Ub_tm1=lambda t :0,␣ (cid:44) → Ub_sp1='1000');H.Project_Boundary_Control();H.Set_Time_Setting(dt);
Now we perform the time integration of the resulting pHDAE system and deduce the behaviour of theenergy variables, the Hamiltonian with respect to the time and space variables, respectively. For thetime discretization, we employ a fully explicit scheme, presented in [45] (Algorithm 2 of Section 4.4) asa first attempt. 28 .Set_Formulation('div')alpha_s, fS, fsig, es, eS, esig, Hamiltonian = H.Integration_DAE();
Figure 6: Hamiltonian (internal energy) versus time for the heat equation.As an illustration, we plot the Hamiltonian function ( i.e. the internal energy) versus time. TheHamiltonian function is constant after seconds, when the boundary control is switched off, as expectedby the first law of thermodynamics. We have provided a general structure for the theoretical and numerical solution of infinite-dimensionalport-Hamiltonian systems. This structure is particularly appealing since PFEM straightforwardly ap-plies. Concerning the numerical solution, PFEM offers the advantage to leverage robust software com-ponents for the discretization of boundary controlled PDEs and time integration.We have applied this strategy on abstract multidimensional linear hyperbolic and parabolic boundarycontrolled systems. We have notably shown that model problems based on the wave equation, Mindlinequation and heat equation fit within this unified theoretical framework. Numerical simulations ofinfinite-dimensional pHs have been performed with the ongoing software project
SCRIMP that has beenbriefly introduced. Finally we have illustrated how to solve three case studies within this framework bycarefully explaining the methodology, and have provided companion interactive
Jupyter notebooks.Beside the generalization of the classes related to the heat and wave equation to the three-dimensionalcase, we plan to propose in
SCRIMP more advanced model problems based on the two-dimensionalShallow Water Equation (SWE) [18, 20], the Kirchhoff model for thin plates [15] and Maxwell’s equations[42]. Furthermore we will investigate both time integration methods that allow structure-preserving timediscretization [32] of finite dimensional pHs and more accurate time integrators for nonlinear pHDAE. Inaddition we plan to enrich the panel of structure-preserving model reduction algorithms to facilitate thesimulation of large-scale port-Hamiltonian systems. This is an essential prerequisite before first attemptsrelated to control design. Further developments foresee the comparisons with well-established algorithmsfor multi-physics problems leading to coupled systems of PDEs.
Acknowledgments
This work is supported by the project ANR-16-CE92-0028, entitled
Interconnected Infinite-Dimensionalsystems for Heterogeneous Media , INFIDHEM, financed by the French National Research Agency29ANR) and the Deutsche Forschungsgemeinschaft (DFG). Further information is available at https://websites.isae-supaero.fr/infidhem/the-project .Moreover the authors would like to thank Michel Salaün and Denis Matignon for the fruitful and in-sightful discussions.
References [1] S. Abhyankar, J. Brown, E. Constantinescu, D. Ghosh, B. Smith, and H. Zhang. PETSc/TS: Amodern scalable ODE/DAE solver library. arXiv:1806.01437, 2018.[2] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring,M. E. Rognes, and G. N. Wells. The FEniCS project version 1.5.
Archive of Numerical Software ,3(100):9–23, 2015.[3] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells. Unified form language: Adomain-specific language for weak formulations of partial differential equations.
ACM Transactionson Mathematical Software , 40(2), 2014.[4] C. Andersson, C. Führer, and J. Åkesson. Assimulo: A unified framework for ODE solvers.
Mathe-matics and Computers in Simulation , 116(0):26–43, 2015.[5] D. Arnold, F. Brezzi, and J. Douglas. Peers: a new mixed finite element for plane elasticity.
JapanJournal of Applied Mathematics , 1(2):347, 1984.[6] D. Arnold and J. Lee. Mixed methods for elastodynamics with weak symmetry.
SIAM Journal onNumerical Analysis , 52(6):2743–2769, 2014.[7] D. Arnold and R. Winther. Mixed finite elements for elasticity.
Numerische Mathematik , 92(3):401–419, 2002.[8] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener,V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes,R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, and H. Zhang.PETSc users manual. Technical Report ANL-95/11 - Revision 3.13, Argonne National Laboratory,2020.[9] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart. Linear port-Hamiltonian descriptor systems.
Mathematics of Control, Signals, and Systems , 30(4):17, 2018.[10] L. Beirão da Veiga, D. Mora, and R. Rodríguez. Numerical analysis of a locking-free mixed finiteelement method for a bending moment formulation of Reissner-Mindlin plate model.
NumericalMethods for Partial Differential Equations , 29(1):40–63, 2013.[11] P. Benner and J. Heiland. Time-dependent Dirichlet conditions in finite element discretizations.
ScienceOpen Research , 2015.[12] D. Boffi, F. Brezzi, and M. Fortin.
Mixed Finite Element Methods and Applications , volume 44 of
Springer Series in Computational Mathematics . Springer-Verlag, Berlin Heidelberg, 2013.[13] A. Brugnoli.
A port-Hamiltonian formulation of flexible structures. Modelling and structure-preserving finite element discretization . PhD thesis, Université de Toulouse, ISAE-SUPAERO,France, 2020.[14] A. Brugnoli, D. Alazard, V. Pommier-Budinger, and D. Matignon. Port-Hamiltonian formulationand symplectic discretization of plate models Part I: Mindlin model for thick plates.
Applied Math-ematical Modelling , 75:940–960, 2019.[15] A. Brugnoli, D. Alazard, V. Pommier-Budinger, and D. Matignon. Port-Hamiltonian formulationand symplectic discretization of plate models Part II: Kirchhoff model for thin plates.
AppliedMathematical Modelling , 75:961–981, 2019.[16] A. Brugnoli, F. L. Cardoso-Ribeiro, G. Haine, and P. Kotyzca. Partitioned finite element methodfor power-preserving structured discretization with mixed boundary conditions. In
Proceedings ofthe 21st IFAC World Congress , volume 53, pages 7647–7652, 2020. Invited session.3017] A. Brugnoli, G. Haine, A. Serhani, and X. Vasseur. Supplementary material for "Numerical ap-proximation of port-Hamiltonian systems for hyperbolic or parabolic PDEs with boundary control". https://doi.org/10.5281/zenodo.3938600 , 2020. Dataset on Zenodo.[18] F. Cardoso-Ribeiro, A. Brugnoli, D. Matignon, and L. Lefèvre. Port-Hamiltonian modeling, dis-cretization and feedback control of a circular water tank. In , pages 6881–6886, Nice, France, 2019. IEEE. Invited session.[19] F. L. Cardoso-Ribeiro, D. Matignon, and L. Lefèvre. A structure-preserving partitioned finiteelement method for the 2d wave equation.
IFAC-PapersOnLine , 51(3):119–124, 2018. 6th IFACWorkshop on Lagrangian and Hamiltonian Methods for Nonlinear Control LHMNC 2018.[20] F. L. Cardoso-Ribeiro, D. Matignon, and L. Lefèvre. A Partitioned Finite-Element Method (PFEM)for power-preserving discretization of open systems of conservation laws. arXiv:1906.05965, 2019.[21] S. Chaturantabut, C. Beattie, and S. Gugercin. Structure-preserving model reduction for nonlinearport-Hamiltonian systems.
SIAM Journal on Scientific Computing , 38(5):B837–B865, 2016.[22] G. Cohen and P. Grob. Mixed higher order spectral finite elements for Reissner–Mindlin equations.
SIAM Journal on Scientific Computing , 29(3):986–1005, 2007.[23] T. J. Courant. Dirac manifolds.
Transactions of the American Mathematical Society , 319(2):631–661, 1990.[24] H. Egger, T. Kugler, B. Liljegren-Sailer, N. Marheineke, and V. Mehrmann. On structure-preservingmodel reduction for damped wave propagation in transport networks.
SIAM Journal on ScientificComputing , 40(1):A331–A365, 2018.[25] G. Golo, V. Talasila, A. J. van der Schaft, and B. Maschke. Hamiltonian discretization of boundarycontrol systems.
Automatica , 40(5):757–771, May 2004.[26] G. Haine, D. Matignon, and A. Serhani. Numerical analysis of a structure-preserving space-discretization for an anisotropic and heterogeneous boundary controlled N -dimensional wave equa-tion as port-Hamiltonian system. arXiv:2006.15032, 2020.[27] A. C. Hindmarsh, P. Brown, K. E. Grant, S. Lee, R. Serban, D. Shumaker, and C. Woodward.SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions onMathematical Software (TOMS) , 31(3):363–396, 2005.[28] P. Joly. Variational Methods for Time-Dependent Wave Propagation Problems. In Mark Ainsworth,Penny Davies, Dugald Duncan, Bryan Rynne, and Paul Martin, editors,
Topics in ComputationalWave Propagation: Direct and Inverse Problems , volume 31 of
Lecture Notes in ComputationalScience and Engineering , pages 201–264. Springer, Berlin, Heidelberg, 2003.[29] R. C. Kirby and T. T. Kieu. Symplectic-mixed finite element approximation of linear acoustic waveequations.
Numerische Mathematik , 130(2):257–291, Jun 2015.[30] R. C. Kirby and A. Logg. Efficient compilation of a class of variational forms.
ACM Transactionson Mathematical Software , 33(3):17–es, 2007.[31] P. Kotyczka.
Numerical Methods for Distributed Parameter Port-Hamiltonian Systems . TUM Uni-versity Press, Munich, 2019. Habilitation.[32] P. Kotyczka and L. Lefèvre. Discrete-time port-Hamiltonian systems: A definition based on sym-plectic integration.
Systems and Control Letters , 133(November):104530, 2018.[33] P. Kotyczka, B. Maschke, and L. Lefèvre. Weak form of Stokes-Dirac structures and geometricdiscretization of port-Hamiltonian systems.
Journal of Computational Physics , 361:442–476, 2018.[34] R. Krug, V. Mehrmann, and M. Schmidt. Nonlinear optimization of district heating networks.
Optimization and Engineering , Sep 2020.[35] M. Kurula and H. Zwart. Linear wave systems on n-D spatial domains.
International Journal ofControl , 88(5):1063–1077, 2015. 3136] Y. Le Gorrec and D. Matignon. Coupling between hyperbolic and diffusive systems: A port-Hamiltonian formulation.
European Journal of Control , 19(6):505–512, 2013.[37] S. Linge and H. P. Langtangen.
Programming for Computations - Python . Springer, 2020.[38] A. Logg, K. A. Mardal, G. N. Wells, et al.
Automated Solution of Differential Equations by theFinite Element Method . Springer, 2012.[39] A. Logg and G. N. Wells. DOLFIN: Automated finite element computing.
ACM Transactions onMathematical Software , 37(2), 2010.[40] V. Mehrmann and R. Morandin. Structure-preserving discretization for port-Hamiltonian descriptorsystems. In , pages 6863–6868, 2019.[41] R. Moulla, L. Lefèvre, and B. Maschke. Pseudo-spectral methods for the spatial symplectic reductionof open systems of conservation laws.
Journal of Computational Physics , 231(4):1272–1292, 2012.[42] G. Payen, D. Matignon, and G. Haine. Modelling and structure-preserving discretization ofMaxwell’s equations as port-Hamiltonian system. In
Proceedings of the 21st IFAC World Congress ,volume 53, pages 7671–7676, 2020. Invited session.[43] R. Rashad, F. Califano, A.J. van der Schaft, and S. Stramigioli. Twenty years of distributed port-Hamiltonian systems: a literature review.
IMA Journal of Mathematical Control and Information ,07 2020.[44] M. Renardy and R. C. Rogers.
An Introduction to Partial Differential Equations . Number 13 inTexts in Applied Mathematics. Springer-Verlag New York, 2nd edition, 2004.[45] A. Serhani.
Systèmes couplés d’EDPs, vus comme des systèmes Hamiltoniens à ports avec dissipation: Analyse théorique et simulation numérique . PhD thesis, Université de Toulouse, ISAE-SUPAERO,France, 2020.[46] A. Serhani, G. Haine, and D. Matignon. Anisotropic heterogeneous n -D heat equation with boundarycontrol and observation: I. Modeling as port-Hamiltonian system. IFAC-PapersOnLine , 52(7):51–56,2019. 3rd IFAC Workshop on Thermodynamic Foundations for a Mathematical Systems (TFMST).[47] A. Serhani, G. Haine, and D. Matignon. Anisotropic heterogeneous n -D heat equation with boundarycontrol and observation: II. Structure-preserving discretization. IFAC-PapersOnLine , 52(7):57–62,2019. 3rd IFAC Workshop on Thermodynamic Foundations for a Mathematical Systems (TFMST).[48] A. Serhani, D. Matignon, and G. Haine. A Partitioned Finite Element Method for the Structure-Preserving Discretization of Damped Infinite-Dimensional Port-Hamiltonian Systems with BoundaryControl. In Nielsen, Frank and Barbaresco, Frédéric, editors,
Geometric Science of Information ,volume 11712 of
Lecture Notes in Computer Science , pages 549–558. Springer, Cham, 2019.[49] A. Serhani, D. Matignon, and G. Haine. Partitioned Finite Element Method for port-Hamiltoniansystems with Boundary Damping: Anisotropic Heterogeneous 2-D wave equations.
IFAC-PapersOnLine , 52(2):96–101, 2019. 3rd IFAC Workshop on Control of Systems Governed by PartialDifferential Equations (CPDE). Invited session.[50] S. Timoshenko and S. Woinowsky-Krieger.
Theory of plates and shells . Engineering societies mono-graphs. McGraw-Hill, 1959.[51] J. Toledo, Y. Wu, H. Ramírez, and Y. Le Gorrec. Observer-based boundary control of distributedport-Hamiltonian systems.
Automatica , 120:109130, 2020.[52] V. Trenchant, H. Ramírez, Y. Le Gorrec, and P. Kotyczka. Finite differences on staggered gridspreserving the port-Hamiltonian structure with application to an acoustic duct.
Journal of Compu-tational Physics , 373:673–697, 2018.[53] M. Tucsnak and G. Weiss.
Observation and control for operator semigroups . Birkhäuser AdvancedTexts: Basler Lehrbücher. Birkhäuser Verlag, Basel, 2009.[54] A. van der Schaft and D. Jeltsema. Port-Hamiltonian Systems Theory: An Introductory Overview.
Foundations and Trends ® in Systems and Control , 1(2–3):173–378, 2014.3255] A. J. van der Schaft and B. Maschke. Hamiltonian formulation of distributed-parameter systemswith boundary energy flow. Journal of Geometry and Physics , 42(1):166–194, 2002.[56] A. J. van der Schaft and B. Maschke. Hamiltonian formulation of distributed-parameter systemswith boundary energy flow.
Journal of Geometry and Physics , 42(1–2):166–194, 2002.[57] A. J. van der Schaft and B. Maschke. Geometry of Thermodynamic Processes.
Entropy , 20(12):1–23,2018.[58] N. M. T. Vu, L. Lefèvre, and B. Maschke. A structured control model for the thermo-magneto-hydrodynamics of plasmas in tokamaks.