Discretization of continuum physics – a comparison of numerical methods from a physical point of view

For numerical calculations in continuum physics partial differential equations and the space-time are discretized. This can be done in different ways. Common approaches are finite difference methods and finite element methods, more rarely finite volume methods are used. Each method has different mathematical properties, which have been discussed in the literature, but they also imply a different physical meaning. This issue is discussed in this article and the connection of finite volume methods to thermodynamics of discrete systems is shown.


INTRODUCTION
There exist many different numerical methods for partial differential equations.Some of these methods can only be applied to special types of differential equations, e.g.elliptical partial differential equations, others are difficult to use in complicated geometries.Numerical methods for partial differential equations use different discretizations of space and time and have different mathematical properties, like convergence and accuracy.These properties have been discussed in the literature.Yet, one problem that is not well discussed is that the choice of the numerical method also implies a certain view about the physical system under consideration, i.e. the physical description is changed from a continuous description to a discrete one.In the following sections finite difference, finite element, and finite volume methods are applied to a model problem and the change in the point of view of physics is discussed.This includes the change of the former continuous space into a discrete one, and the continuity of the solution and also of the (numerical) fluxes.

A model problem
As numerical calculations for the systems of partial differential equations, i.e. the balance equations, are very time-consuming and difficult to handle, we will consider a simple model problem for exemplary calculations: the Laplace Problem.
A numerical solution U is calculated, which approximates the (exact) solution u to the two-dimensional Laplace Problem for mixed boundary conditions: where Ω ⊂ R 2 is a bounded Lipschitz domain with polygonal boundary Γ.Dirichlet conditions are assumed on a closed subset Γ D of Γ, while Neumann boundary conditions are assumed on Inhomogeneous Dirichlet conditions are incorporated by using the

FINITE DIFFERENCE METHODS
Finite difference methods (FDMs) [1] emphasize the grid/mesh points.Using FDMs one calculates a solution on a regular mesh; the approximated solution exists only pointwise.The construction of a continuous solution is a different matter, which can be tried afterwards either by fitting by a global function or interpolation (for example, by using splines).

The discretization of the model problem
The discretization with FDMs is based on the simplest approach and probably the first one that comes to mind: replacing the differential operators by difference quotients.
With the help of symmetric differences the discrete Laplace operator can be constructed in the following way.First write the Laplace operator as a sum of second order Cartesian derivatives: then discretize each second derivative by using symmetric differences with x j−1 = x − h and x j+1 = x + h and sum up: The result is a so-called five-point-star: The result of an FDM is a set of approximate values of the function on a discrete set of points.Often the result is displayed with interpolating lines, which is wrong!The interpolation is a different procedure for which new assumptions have to be made; one has to be careful that these assumptions are not in contradiction with assumptions made for applying the FDM.

FINITE ELEMENT METHODS
Finite element methods (FEMs) [1] also emphasize the grid-nodes, but in addition take care of the interpolation.A remarkable point is that the system of partial differential equations is not solved directly, but instead its weak formulation is solved.There exist different types of FEMs according to the type of the mesh and the ansatz-functions used: • P 1 : triangular mesh, linear ansatz-functions; • P 2 : triangular mesh, quadratic ansatz-functions; • Q 1 : quadrilateral mesh, (bi-)linear ansatz-functions; • Q 2 : quadrilateral mesh, (bi-)quadratic ansatz-functions; • . . ., including mixed FEMs.
The weak formulation of this problem reads: 3.1.The discretization of the model problem: Galerkin discretization for a The model problem is discretized using a standard Galerkin method for implementation [2], where H 1 (Ω) and H 1 D (Ω) are replaced by subspaces S and Let (η 1 , . . ., η N ) be a basis of S and (η i 1 , . . ., η i M ) be a basis of S D , where Then the discretized problem (10) is equivalent to Furthermore, let which means we choose the same ansatz for the solution as for the test function.Then the discretized problem results in a linear system of equations Ax = b, where the coefficient matrix and the right-hand side are defined as The coefficient matrix is symmetric and positive definite.So the linear system of equations has exactly one solution x ∈ R M , which gives the Galerkin solution  We assume a polygonal boundary Γ of the domain Ω; then Ω can be covered by a regular triangulation T of triangles and quadrilaterals, i.e.Ω = ∪ T ∈T T with each T being either a closed triangle or closed quadrilateral.The nodes N of the mesh lie on the vertices of the quadrilaterals or triangles, whereas no node lies on an edge of an element.The elements do not overlap and each edge E ⊂ Γ of an element T ∈ T belongs either to ΓN or to ΓD .
The coordinates of the vertices of an element determine the local stiffness matrix.For a triangular element T with the vertices (x 1 , y 1 ), (x 2 , y 2 ), and (x 3 , y 3 ) and the corresponding basis functions η 1 , η 2 , and η 3 , holds true.See Figs 1 and 2 for ansatz-functions over a triangle and quadrilateral and Fig. 3 for a patch of ansatz-functions.

TAKING ONE STEP BACK: HOW THE LOCAL BALANCES ARE DERIVED/ MOTIVATED?
The first two methods are mathematical approaches to the problem of finding an approximate solution to the system of partial differential equations.A more physical approach can be motivated by remembering that the system of partial differential equations under consideration is actually a set of balance equations and how the local balances have been derived/motivated [3].The global balance equation for an arbitrary extensive quantity with specific density ψ, flux J ψ , production density π ψ , and supply density σ ψ reads: where da = nda, and n is the outward unit normal vector to the surface ∂ G.As an example we consider the balance of mass.Then we have ψ ≡ 1 and σ 1 ≡ 0, π 1 ≡ 0, which expresses the conservation of mass: Applying the Reynolds transport theorem to the time derivative of the integral in equation ( 18) and the Gauss theorem to the flux over the boundary in equation ( 18), we obtain where w is the mapping velocity of the control volume.
As these equations hold true for arbitrary control volumes, one can conclude that the integrand must vanish.This results in local balance equations of the form This is exactly what finite volume methods (FVMs) are based on.

FINITE VOLUME METHODS
Finite volume methods [1] in one sense reverse the derivation of the local balance equation and are therefore the most physical approach of the mentioned methods.There are many different types of FVMs, some focusing on the mesh(-nodes) and interpolation and some focusing on the control volumes.Finite volume methods are used if differential equations have a divergence form, i.e. contain parts like with K : Ω → R d,d , c : Ω → R d , and r, f : Ω → R, a part which is also included in parabolic equations of the form which is a generalized form of balance equations.

The discretization of the model problem
Our stationary model problem reads with k ≡ 1.The starting point is to divide the area into control volumes, which have a polygonal boundary, and integrate the differential equation over the control volumes: Here Λ := i ∈ Λ|a i ∈ Ω , a i is an inner node of the associated Delaunay triangulation of Ω, and Ω i is the Voronoi polygon (Wigner-Seitz cell) of a i (see Figs 4,5).Then, by using the Gauss law and ∂ Ω = ∪ j∈A i Γ i j , we obtain The approximation of the normal derivative which represents the flux is the crucial point.As an example difference quotients are used here, which leads to a finite difference-like scheme where , m i j is the length of Γ i j , and Using as an approximation for the integral over f , we finally arrive at

INTERMEDIATE COMPARISON OF THE METHODS
In previous sections the discretization of the continuum for numerical methods was discussed.Figure 6 shows an illustrative comparison of how different methods "approximate" the solution.
As Figure 6 shows, the results of different numerical methods reveal significant differences.By using a finite difference method, the solution is approximated only on a discrete set of points.Finite element methods calculate the approximate solution at discrete points and the interpolation between the points.Finite volume methods result in an approximate (averaged) value, which belongs to the computational cell as a whole.To summarize: FDM: solution and flux are non-continuous, exist only on a discrete set of points; FEM: solution exists on a continuous space, fluxes show jumps; FVM: solution exists on a continuous space, but has jumps, fluxes are continuous.
The behaviour of FVMs to assign the averaged value to the computational cell and the calculation method by use of fluxes over the boundary of the cell shows a significant similarity to thermodynamics of discrete systems.

THERMODYNAMICS OF DISCRETE SYSTEMS
In thermodynamics of discrete systems the system is described as a whole; no internal structure or distribution is taken into account.Exchange with the environment is due to fluxes over the boundary (see Figs 7 and 8).
Different situations in thermodynamics of discrete systems: One equilibrium system in contact with a reservoir: This situation is the standard case discussed in thermodynamics books.It has been well known for a long time.See Fig. 7a.
One non-equilibrium system in contact with a reservoir: This situation is also well understood, but usually not discussed in books, an exception is [4].See Fig. 7b.
Two or more non-equilibrium systems in contact with each other and a reservoir: This situation has only been discussed recently in [5,6].It is not mentioned in books at all.A generalization from two systems to more should be possible by partial contact quantities.See

Contact quantities and replacement quantities
The familiar quantities from thermostatics like the (thermostatic) temperature, pressure, and chemical potential do not exist in non-equilibrium situations.However, non-equilibrium counterparts can be defined 1 (see [4]): the contact temperature Θ, dynamic pressure p, and dynamic chemical potential µ.
These non-equilibrium contact quantities are defined by the exchange of the system with its equilibrium environment by defining the inequalities For a discussion of these inequalities see [4].
The situation changes if a contact between two non-equilibrium systems is considered.The first idea that one might have is to use the contact temperature of the non-equilibrium environment, but this does not work.Instead, the non-equilibrium environment is replaced by an equilibrium environment which causes the same fluxes.The thermostatic temperature of the replaced environment is called the replacement temperature θ * , and analogue for the pressure π * and chemical potential ν * .The defining inequalities read: The difference is that in this case the replacement quantities (θ * , π * , and ν * ) and not the contact quantities are defined by the inequalities.

Compound deficiency and excess quantities
A system can be composed of subsystems.Of course, the accuracy depends on how the system is described -as a whole, without taking the composition into account or as a compound system.This difference is the compound deficiency; it can be taken into account by introducing excess quantities [5,6].
At the boundary of a subsystem the actual value of the quantity is the averaged quantity plus the excess quantity (Fig. 9).
In this way excess quantities can be used to formulate fluxes and jump conditions between computational cells. 1 At least as long as one non-equilibrium system in contact with a reservoir is considered.

CONNECTION BETWEEN FINITE VOLUME METHODS AND THERMODYNAMICS OF DISCRETE SYSTEMS
Each of the control volumes of the FVM forms a discrete system (Schottky system) in the sense of thermodynamics.The value of the variable belongs to the system as a whole and the exchange of quantities is described by fluxes over the boundaries of the system.One point where a good knowledge of thermodynamics of discrete systems can be useful is the approximation of the flux over the boundary An example of the usefulness of thermodynamical considerations in computational algorithms is the derivation of a second order accurate algorithm for hyperbolic first order partial differential equations, which is stable for the Courant number ∆t/c∆x ≤ 1 in [7] -it makes use of excess quantities to define jump relations at the boundary between cells.This algorithm is equivalent to the wave-propagation algorithm proposed by Bale et al. in [8], who actually proved the stability and accuracy.

CONCLUSIONS
Different numerical methods for partial differential equations were discussed.A correspondence between FVMs and thermodynamics of discrete systems was shown by introducing interacting non-equilibrium subsystems.The value of variables belongs to the discrete system (control volume) as a whole.This leads to jumps at the interfaces between control volumes and a stepwise approximation of the solution using an FVM.Compared to the shape functions used for an FEM, this approximation seems rude, but the advantage is that one can handle the fluxes between the control volumes.Finite element methods on the other hand give a steady approximation of the solution, but the value of the variables varies across the cells; therefore a correspondence to discrete systems is not given.

Fig. 6 .
Fig. 6. Results of different methods.(a) An FDM result: Approximative solution at points.The continuum has been eliminated.(b) An FEM result: Approximative solution at points with interpolation.The continuum has been preserved, no jumps are allowed.(c) An FVM result: Averaged value over the control volume.The continuum has been preserved, jumps at interfaces between control volumes are possible.

Fig. 8 .Fig. 7 .
Fig. 7. Two discrete systems in contact with each other and with the environment.(a) Equilibrium.(b) Non-equilibrium.

Fig. 8 .
Fig. 8. Two discrete systems in contact with each other and with the environment.