Systematisation of Systems Solving Physics Boundary Value Problems

A general conservation law that deﬁnes a class of physics ﬁeld theories is constructed. First, the notion of a general ﬁeld is introduced as a formal sum diﬀerential forms on a Minkowski manifold. By the action principle the conservation law is deﬁned for such a general ﬁeld. By construction, particular ﬁeld notions of physics, e.g., magnetic ﬂux, electric ﬁeld strength, stress, strain etc. become instances of the general ﬁeld. Hence, the diﬀerential equations that constitute physics ﬁeld theories become also instances of the general conservation law. The general ﬁeld and the general conservation law together correspond to a large class of relativistic hyperbolic physics ﬁeld models. The parabolic and elliptic models can thereafter be derived by adding constraints. The approach creates solid foundations for developing software systems for scientiﬁc computing; the unifying structure shared by the class of ﬁeld models makes it possible to implement software systems which are not restricted to certain predeﬁned problems. The versatility of the proposed approach is demonstrated by numerical experiments with moving and deforming domains.


Introduction
In this paper we focus on second-order boundary value problems (BVP's) related to physical field theories.BVP's and their numerical solution methods is an extensively studied field of science.Still, many practical challenges remain, e.g.: i) One may have a problem to which there is no software system available.ii) The software systems are laborious if not hard to extend beyond their original purpose and such extensions increase the complexity of the system.iii) In case of incorrect results, it is tedious to distinguish between simple user errors and errors in reasoning.iv) Users often have to learn many software specific details.

Differential geometric models
In the late 1990's and early 2000's, Bossavit et al. developed and introduced the so called "geometric approach" into electromagnetism [7].In addition, in 1997 the idea of a "discrete Hodge operator" [39] was introduced to reveal the key mathematical structures behind finite difference and finite element kind of methods [9].At the time the finite difference method [37] was commonly explained in a rather elementary manner in Cartesian coordinate systems following K. Yee's original paper [43] from 1966.Bit later the scientific community in elasticity picked the idea and the geometric approach became known also as "Discrete exterior calculus" (DEC) after Hirani [18].
We have further developed the geometric approach and created a generic software system based on it.The system can be employed to solve hyperbolic application problems from classical and quantum physics [24], [29], [31], such as electromagnetic, elastic, and acoustic wave problems, the Schrödinger equation [13], or Gross-Pitayevskii equations [30], and so on.We explain the mathematical foundations of the software system in [22].The implementation of the simulation software and the various mesh structures which we have employed are described in detail in [29].
To explain the methodology, we will first outline a theory of ordinary gauge theories on form bundles. Thereafter we will briefly discuss the extension to Clifford and tensor algebra.Exterior (or Grassmann) algebra [14], [23] is the Clifford algebra [17], where the quadratic form is identified to zero, and Clifford algebra itself is a quotient algebra of tensor algebra [26].We assume a Minkowski manifold [3], [15], and describe the proposed methodology in steps from the foundations.

Formal sums of field configurations
The field-notion in physics involves an idea of assigning numbers to geometrical objects of space-time, such as to points, (virtually) small segments of oriented lines, etc.These numbers represent observations made by measurements, and they can be interpreted as the values differential forms yield on p-vectors [5].
Let us start from ordinary differential forms, which come with a degree from zero to the dimension n of the manifold.Since we are not after any particular field configuration, forms of a particular degree are not in our interest.We hide the information of the degree by introducing a formal sum of differential forms of all degrees: where α p ∈ {0, 1}.Please notice that with ordinary differential forms and in the n-dimensional case the number of p-forms is n p and the formal sum has the total of 2 n degrees of freedom.
By operating with F the emphasis is shifted from particular degrees to the property that all forms map some p-vector, 0 ≤ p ≤ n, to scalars.

Differentiation and the action principle
Next, we need to introduce differentiation for F. This is straightforward as smooth p-forms are differentiated with the exterior derivative d, and so is also F. Ordinary gauge theories are characterized by pairs of differential equations, such as electromagnetic theory [36] is described by Maxwell's equations.The gauge-theoretic view is that differential equations follow from the action principle [3].An action is the integral of a Lagrangian L over a manifold, and differential equations correspond to the critical points of the action.
A large class of models in ordinary gauge theories have to do with the conservation of some quadratic notion.We equip the Minkowski manifold with a metric tensor providing us also with a Hodge operator .Then, we assume F is an exact field, F = dH where H = h 0 + . . .+ h n is a potential.In addition, for the source terms we introduce another formal sum G = g 0 + . . .+ g n .Now, an action of the desired type can be given by

The differential equations are then obtained as follows. The variation of action
where H α = H + αδH and by insisting on the variation δ A to vanish for all δH yields the critical points of A and the corresponding differential equations.Hence, the action principle implies that at all (ordinary) points on the Minkowski manifold the following differential equations dF = 0 and d F = G should hold.These equations can be expressed as the diagram in Figure 1.

Fig. 1 Diagram of differential equations
Let us also express the action principle as a diagram.For brevity, to introduce such a diagram we assume G to vanish.Then, the Lagrangian of the action becomes L = n p=0 1 2 f p ∧ f p , and the definition of the Hodge operator implies, that each component L p satisfies where ω 0 is the unit n-volume of Minkowski space and q p is the quadratic refinement of the Minkowski bilinear form •, • .The L p 's form a product space L = D f × D * f equipped with projections π f : L → D f and π * f : L → D * f satisfying the following universal property: For every action A and Lagrangian L there is a unique map a : A → L and l : L → L such that the diagram of Figure 2 is commutative.

Fig. 3 Diagram DGOrd
The combination of the two diagrams of Figure 1 and Figure 2 results in a diagram presenting how the action with the Lagrangian defines differential equations for a pair of fields, which are in a Hodge relation to each other.We call this diagram by the name DGOrd (designating that it involves ordinary differential forms), simplify it a bit -object L is left out-and draw it in Figure 3. Ordinary gauge theories include also other type of differential forms than ordinary ones, which are also essential in mathematical physics.For instance, in elasticity [1] E-valued forms, vector and covector-valued forms [3], [16], [23] are needed [21], [32], [33] [34], [35].
Let us next extend the idea of formal sums of differential forms to E-valued forms.By construction, such formal sums of E-valued forms can be differentiated with the exterior covariant derivative d ∇ , where ∇ is the connection.To introduce the Lagrangian as a quadratic refinement of the Minkowski metric, the Hodge operator should be extended to E-valued forms such that L = n p=0 1 2 f p ∧ f p becomes a formal sum of scalars.We denote such a Hodge operator by E .In the same manner End(E)-valued [3] (i.e., matrix valued) forms are needed for example in Yang-Mills theory [2], [38], [41], [42].The Hodge operator End is now extended to End(E)-valued forms so that the Lagrangian becomes a formal sum of scalars.
Formally, there exists an abstract diagram DGA shown in Figure 4, and mappings M 0 : DGA → DGOrd, M E : DGA → DGE and M End : DGA → DGEnd.They map the abstract diagram DGA to more concrete diagrams DGOrd of ordinary forms, DGE of vector valued forms and DGEnd of matrix valued forms.They also map the hodge duality to operators , E , and End , respectively.M O maps the differentiation to exterior derivative d.M E and M End map it to d ∇ .
This construction suggests that mappings M O , M E and M End represent various models of the "theory of differential geometric models" represented by DGA.This is a step towards a category theoretic representation of physics field theories.
Remark: Hyperbolic wave problems in physics are particular examples of our models.DGA can also be concretized to elliptic and parabolic models.Later, as an example, we show how to concretize the Schrödinger equation from the general setting.As we use differential geometric formalism, the canonical way to discretize all the considered models is to use DEC.

Concretization of particular models
Next, to verify the usefulness of the theory and its models in scientific computing, let us exemplify how particular models are concretized from the theory.This also highlights the pragmatic significance of a proper mathematical theory; resources become more efficiently exploited, if software systems are designed to realize theories instead of particular models.
Let us start by concretizing DGOrd to four dimensional differentiable manifold Ω with Minkowski metric, signature (−, +, +, +), and a decomposition of space-time into space and time-like components; Ω = Ω t × Ω s .Symbols F and G denote formal sums of p-forms, F = f 0 + . . .+ f 4 and G = g 0 + . . .+ g 4 and consequently, differential equations dF = 0 and d f = G can be written as: . By i) decomposing p-forms into time-like components and only space-like components, ii) and exterior derivative to space and time-like components, and iii) applying the Leibniz rule, we obtain an equivalent system of Equation ( 1) [22].
Here d and are now the exterior derivative and the Hodge operator, respectively, in the space-like component Ω s of manifold Ω. Subscript s in the f p and g p 's denotes the space-like component f p s of (p + 1)-form f t dt ∧ f p s .This system of equations and its natural transformations cover a wide class of physics field theories.By construction, all the models covered by the theory are relativistic.Each particular model corresponds to a choice of F and G as demonstrated next with some examples.
For Maxwell's equations [36] in space and time, F is chosen to be the Faraday field and G the source charges q and currents j [3]: and G = j − dt ∧ q.This corresponds to setting f 1 s = −e, f 2 = b, g 1 = j, and g 0 s = − q, and by substituting these to the system of Equation (1), we obtain We have considered permittivity and permeability µ as properties of Ω s .Thus, they are embedded into the Hodge operators and µ [6].The non-relativistic Schrödinger equation [13] can also be concretized from DGOrd by imposing some simplifying constraints on the general model.For this, we first set − ϕ R = ϕ I and choose where is the reduced Planck constant and m is particle's mass.By substitution to the general system one obtains: The relativistic property is next lost by a modelling decision.The term ∂ t q R is assumed to vanish.This implies that also −∂ t q I = ∂ t q R = 0. Now the middle equations become tautologies and the system is reduced to the pair By mapping differential forms to vector fields and mapping the exterior derivative to the corresponding differential operators of vector analysis, the textbook version div grad ϕ = −iV ϕ results.It is defined using complex arithmetic which restricts it to flat Minkowski manifold only.The relativistic intermediate stage obtained from the general model can also be implemented on curved space-time.This is an interesting topic to be numerically tested.
Small-strain elasticity is naturally modelled using E-valued forms.Recall that in this case the differential equations on the manifold Ω t × Ω s take the form d ∇ F = 0 and E d ∇ E F = E E G. Analogously to previous, Leibniz rule and the spacetime split of forms and exterior covariant derivative d ∇ results in structurally similar general system as in the case of ordinary forms.The exterior derivative d is simply replaced by d ∇ and is replaced by E .
The model of elasticity now arises by the choice f 0 s = u, f 1 = ε, g 0 = − f v , where the vector-valued 0-form u is the time-derivative of displacement ν, u = ∂ t ν, the vector-valued 1-form is linearized strain, and g 0 is the source force term.Substituting this choice back to the system of equations yields The Hodge operator C contains the parameters of the stress-strain relation, and density ρ is embedded to ρ .Since u = ∂ t ν, the first equation states that = d ∇ ν and the second equation is automatically satisfied.As a result, we get the elasticity equations which are, for convenience, written out also in Euclidian space and using vector analysis notation: The final example is Yang-Mills equations [42] where the field configurations are End(E)-valued forms.As Yang and Mills developed their theory as an extension to Maxwell's theory, Yang-Mills equations are concretized from the system of differential equations for End(E)-valued forms in the same manner as Maxwell's equations are concretized from the system for ordinary forms.Such a process results in Clifford and tensor algebra.The "theory of differential geometric models" presented is not complete for the needs of software design: First, tensor algebra is the most general algebra for vector spaces over scalars, and all field configurations share the structure of a vector space.Second, a Clifford algebra is unital associative algebra generated by a vector space equipped with a quadratic form [17]. Third, exterior algebra is the Clifford algebra when the quadratic form is zero.Clifford algebra seems to provide us with a better starting point as certain Clifford algebras, such as Pauli or Dirac algebra [15], are very important in mathematical physics.
We construct the universal Clifford algebra as a subalgebra of the algebra of linear transformations [17].Let F be a scalar field and denote In the algebra L(Λ(V )) of linear transformations of Λ(V ) map M v is defined as the linear extension of M v (1) = v, and where v k denotes the term to be omitted from the product, and where B(•, •) is the Minkowski bilinear form.Thus, M v is exterior multiplication by v and δ v is interior multiplication with respect to the inner product induced on Let the (metric compatible) covariant derivative be mapped by functor C from the tensor bundle to the Clifford bundle.The image of the covariant derivative in the Clifford bundle is denoted by ∇.The codomain of map ∇ : L(Λ(V )) → L(Λ(V )) can be decomposed into components corresponding to the exterior and interior multiplication, and consequently we may write cod(∇) = cod(∇ e ) ⊕ cod(∇ i ).If the covariant derivative is mapped to the exterior bundle with functor D, then the image of the covariant derivative is d ± d , where the sign depends on grade p and dimension n.
Tensor bundles and tensor algebra provide us with a starting point general enough for the theory needed in software design for ordinary gauge theories.The theory should not, however, be tied to the category of sets.We seek for a category which just condenses the essentials of differentiation, of the metric properties of spacetime, and of the action principle.Software based on such a theory is not bound to any specific algebra.This enables the end users to employ algebras that fit best with their needs.

Some numerical experiments with space-time models
In our earlier papers, we have demonstrated the proposed approach with several numerical experiments.Such experiments include simulations with acoustic, elastodymanic, electromagnetic and quantum mechanic waves [28], [29], [30], [31].For the extended accuracy, the mesh structures play an essential role [29].The numerical scheme can also be optimized by locally adaptive time-stepping and by tuning the dicerete Hodge operator, e.g., for time-harmonic waves.In certain cases such optimizations can improve the efficiency of the simulation even by orders of magnitude as reported in [28].
The formulation of the general model in Minkowski space provides additional benefits.It is namely possible to simulate the wave propagation in moving (and even deforming) spatial domains.In the papers [28], [29], [30] and [31], the spatial mesh generation together with the associated spatial finite difference approximation and time-stepping were considered as separate entities, without emphasizing the fact that the usual leap-frog time integration scheme for first order systems could also be derived from geometrical principles analogous to the spatial mesh generation which is based on the Delaunay-Voronoi duality.This chapter contains numerical experiments that demonstrate how the general model of the Minkowski space can be discretized.The construction of the spacetime model begins by creating a mesh that fills the entire space-time domain of the computation.When generating a mesh, one should ensure that a valid dual mesh is available.The dual mesh is made up of cells each having an orthogonal counterpart in the (primal) mesh.Orthogonality is defined such that any vector from a primal cell has zero inner product (in the sense of Minkowski metric) with any vector from the corresponding dual cell.

Transforming computational cavity
Figure 5 illustrates the solution of the one-dimensional time-dependent wave problem in a moving cavity.We build a simplicial mesh in a two-dimensional space with one spatial axis and a time-axis.Then we attach a floating point number to each primal 1-cell (edge) to construct a discrete version of F including only 1-form term.The initial values are set at time t = 0 (at the bottom of the figure) to trigger a wave pulse.Elsewhere in the computing domain, the values of F are explicitly solved by following the equation dF = 0. Since the dual mesh and the discrete Hodge operators are constructed using Minkowskian metric, the solution is a traveling wave with a propagation speed of 1 in both directions (see the right-hand-side of the Figure 5).Figure 6 shows a numerical experiment where the same approach has been applied to solve a two-dimensional time-dependent wave problem in a rotating cavity.In this case, the mesh is three-dimensional and the shape of the two-dimensional base mesh (spatial cross-section of the space-time mesh) resembles a boomerang.The space-time mesh is twisted around the time axis, causing the cross-section to rotate as time progresses.The field F to be solved is a 2-form which is discretized by attaching one floating-point number to each 2-cell (face) of the mesh.By initializing F as an impulse at the initial time and solving dF = 0 inside the computational domain, we detect a wavefront propagating at speed 1 and reflecting from the moving walls.A video of this numerical experiment can be found at the following url: https://urly.fi/1oxH.To prove the generalizability of the method, we present yet another experiment where we solve a three-dimensional acoustic-like wave problem in a shrinking computational domain.We create a (3+1)-dimensional simplicial mesh that, at time t = 0, fills a three-dimensional spatial volume as illustrated on the left of Figure 7.The element lengths of the mesh are proportional to the term 1 − 0.3t.This means that the element sizes decrease exponentially in terms of the number of time steps.The point (0, 0, 0, 1  3 ) T of convergence is never reached in the simulation.In order to reduce the amount of memory required, the mesh duration over time is chosen as short as possible.We integrate 1-form F over mesh by explicitly solving dF = 0.When integration over mesh is completed, the last calculated terms are copied as the initial values of the next iteration and the integration is repeated.In this way, the task can be integrated as long as desired, without having to store the entire mesh in memory.The resulting field of time-integration is illustrated in Figure 7.A video of this numerical experiment can be found online at url: https://urly.fi/1oWx.

Local time-stepping and stability
Traditionally, the Courant-Friedrichs-Lewy (CFL) condition sets an upper limit for the length of maximal time step.The smaller the spatial element size is, the shorter the time step must be in order to achieve numerical stability.When the spatial element length is not constant, local time-stepping can speed up the integration of time-dependent wave problems.This section shows how to create local time-stepping methods using the space-time integration.
Let's start with a (1+1)-dimensional example and consider a one-dimensional spatial mesh consisting of unevenly distributed nodes and line-segments (edges) between them.Nodes of the spatial mesh are copied at regular intervals in the time direction using individual step sizes ∆t for the nodes.The length of the time step is set to the maximum length that obeys the inequality ∆t < c∆x, where ∆x is the length of the shortest neighboring edge and c is a constant.The space-time structure is completed as the Delaunay mesh.The mesh is 2 units wide in spatial direction and 1 unit high in time-direction.The mesh and its dual mesh are illustrated in Figure 8.The field under consideration is a 1-form and it is formatted and integrated in the same way as in the previous section.When the integration over the mesh is completed, the last calculated terms are copied as the initial values and the integration is repeated.In this way, we are able to reuse the mesh again and integrate over time as long as necessary.
We consider stability of the time-integration in long term simulations with two different constants c = 1.0 and c = 0.9.The results are illustrated in Figure 9.The conclusion is that the time integration is not stable with the constant of c = 1.0.The noise in the resulting field is visible already after 50 iterations.However, with the constant of c = 0.9, the system is stable because no dispersion is detectable even after 200 iterations.The condition used for the time step length seems to be a good first guess to replace the CFL condition in the asynchronous space-time integration.
We also investigate local time-stepping in a (2+1)-dimensional wave problem.The mesh is constructed by creating a two-dimensional circular base mesh with varying element sizes.We limit the individual time step ∆t of each node by the relation ∆t < c∆x, where ∆x is length of the shortest spatial edge next to the node.The structure of the space-time mesh is determined as a Delaunay mesh and a truncated mesh is visualized in Figure 10.
We integrate 1-form over time and consider the numerical stability of long-term simulations.From Figure 11, we find that by limiting the time step length with the  constant of c = 0.7, numerical stability is not achieved.We observe noice in the resulting field already after 5 units of time.Instead, using the constant c = 0.6, we keep the integration stable and do not observe any dispersion in the resulting field even after 100 units of time.

Conclusions
In this paper we have considered the common structure of boundary value problems.
The structure is based on ordinary gauge theories on form bundles.We have presented models from classical and modern physics as particular examples of the system.The finite dimensional models are constructed with discrete exterior calculus (DEC) method from the models expressed with differential forms.The computational mesh is based on Delaunay-Voronoi duality with Minkowskian metric.With the corresponding construction of discrete Hodge operators, it enables numerical simulations in moving and deforming spatial domains.Adaptive time stepping can also be implemented by utilizing the geometry of space-time mesh.Numerical results show that the software system based on the systematic structure is applicable in boundary value problems in one, two, and three spatial dimensions.

Fig. 2
Fig. 2 Diagram of the action principle

Fig. 5 A
Fig. 5 A space-time approach to solve of a time-dependent wave problem in a moving cavity: The mesh with simplicial cells (purple edges) and corresponding Minkowski dual mesh (blue edges) are illustrated on the left.The solution of a wave problem is shown on the right.The color components red and green correspond to dx and dt components of the resulting 1-form.The colors are normalized such that grey indicates the zero field.

Fig. 6
Fig. 6 Simulation of wave propagation in rotating 2-dimensional cavity: The (2+1)-dimensional space-time mesh is illustrated on the left.The red color (dark) at the bottom indicates the past time and the cross section of the mesh at the current time is shown on the right.The color components red, green, and blue represent the components dx ∧ dy, dx ∧ dt, and dy ∧ dt of the resulting 2-form, respectively.The figure is normalized such that the grey color indicates the zero field.

Fig. 7
Fig. 7 Wave propagation in a shrinking 3-dimensional cavity: Cross-section of the space-time mesh and dt-component of the resulting field are presented at five instances of time.

Fig. 8 AFig. 9
Fig. 8 A (1+1)-dimensional mesh with variable spatial edge lengths ∆x.The condition for time step size is ∆t < 1.0∆x.Primal and dual edges are illustrated with purple and blue colors, respectively.

Fig. 10 A
Fig. 10 A (2+1)-dimensional mesh with variable spatial edge lengths ∆x and with condition ∆t < 0.6∆x for the time step size.

Fig. 11
Fig. 11 Cross-sections of fields at different time instances and under different conditions for a time step length.The color components red, green, and blue correspond to dx, dy, and dt components of the resulting field, respectively.The grey color indicates the zero field.