On the Verification of the Pedestrian Evacuation Model

In this article we deal with numerical solution of macroscopic models of pedestrian movement. From a macroscopic point of view, pedestrian movement can be described by a system of first order hyperbolic equations similar to 2D compressible inviscid flow. For the Pedestrian Flow Equations (PFEs) the density ρ and the velocity v are considered as the unknown variables. In PFEs, the social force is also taken into account, which replaces the outer volume force term used in the fluid flow formulation, e.g., the pedestrian movement is influenced by the proximity of other pedestrians. To be concrete, the desired direction μ of the pedestrian movement is density dependent and is incorporated in the source term. The system of fluid dynamics equations is thus coupled with the equation for μ. The main message of this paper is the verification of this model. Firstly, we propose two approaches for the source term discretization. Secondly, we propose two splitting schemes for the numerical solution of the coupled system. This leads us to four different numerical methods for the PFEs. The novelty of this work is the comparative study of the numerical solutions, which shows, that all proposed methods are in the good agreement.


Introduction
Modeling of the movement of pedestrians is an interesting non-linear complex problem, where physical and environmental conditions (e.g., visibility see [1]), as well as social and psychological conditions (e.g., fatigue, see [2]), manifest themselves. The complexity of the dynamics of pedestrian movement is illustrated by the emergence of collective effects and self-organisation phenomena. One of the most important practical applications of these models is the ability to simulate the movement of a crowd in panic situations, such as floods, fires, etc. In such cases, it is necessary to have evacuation plans in place and try to predict people's behavior.
In general, there are two classes of the pedestrian behaviour models, macroscopic and microscopic. We are mainly interested in the macroscopic approach where pedestrians are considered as a continuum and their behaviour is described by two-dimensional fluid dynamics equations for an inviscid adiabatic flow subject to outer volume forces. This approach was proposed in [3].
The outer volume forces represent a relaxation term towards a desired velocity. Since the model is devised for room or open space evacuation situations, the desired pedestrian velocity depends on the density distribution. The pedestrians want to avoid the regions with high densities, so the magnitude of the desired velocity is a decreasing function of the density. Such a function stems from the experimental investigations of pedestrian evacuation problems. For various definitions of density dependent magnitude of desired velocity see [4][5][6]. A detailed overview can be found in [7]. Also the desired direction of pedestrian velocity is density dependent. The desired direction of a pedestrian motion is given as a unit tangential vector to the fastest path to the exit. This works like a car navigation system. The pedestrian knows the magnitude of the velocity in which he would like to go and seek the fastest path to exit. Instead of seeking the fastest path, the so called Eikonal equation can be used for the determination of the desired direction of pedestrian velocity. See [8] for the detailed explanation. The above described model falls into the class of macroscopic models, where the state of the system is described by locally averaged quantities. It was originally proposed in [9] as an extension of the vehicular traffic models described in [10,11]. More models originated from different pedestrian behaviour and strategies are proposed and analyzed in [7]. The reader will find a classification of macroscopic models, and some qualitative properties such as hyperbolicity are further studied. The model of pedestrian with disctinct desired direction is studied in [12], here the two-phase pedestrian flow model is taken into account.
In the microscopic approach, the pedestrian is considered as an individual interacting with others pedestrians. Example of this approach is the model, where the pedestrians are described by the system of N ordinary differential equations (ODE), here N stands for the number of pedestrians. This model is studied in [13]. Here are stated simulations of the high-density crowd behaviour in the real world objects and situations such as sporting events and pilgrimage places. The transition from microscopic (ODE) to macroscopic model is listed in [14]. Reader founds deep compilation of the literature devoted to booth approaches in [3]. Here is stated more detailed description and differentiation of microscopic models such as a social force model, cellular automata models (CA) or gradient flow methods (macroscopic approach). The comparison of the cellular automata model with the PFEs (macroscopic) can be found in [15], where the new Cellular Automaton (CA) model is proposed for the simulation of the pedestrian flow and the comparative study is presented. The level between microscopic and macroscopic scale is represented by the mesoscopic scale. An example of this kind of models is kinetic model described in [16].
It is seen from a review of the literature that extensive research is currently underway on modelling pedestrian evacuation and related situations. In this paper, we are interested in a model that is based on fluid dynamics equations and formulated as a hyperbolic system with a source term on the right hand side of the equation. In accordance with [17], such a system is usually solved by the splitting method. Because there is a lack of experimental data from real-life situations in this case, it is difficult to validate that the numerical results correspond to real-world situations. But what is possible is to study whether the numerical solution is correct.
To this end, we have proposed two substantially different numerical methods of splitting: Splitting Scheme I and Splitting scheme II. The paper emphasizes which parts of discretization using the Finite Volume Method (FVM) and the combined Finite Volume and the Discontinuous Galerkin method (FV-DG), are original. The choice of the numerical flux in both methods is original (the original numerical flux of the Vijayasundaram type, the method of solving the Eikonal equation using Dijkstra's algorithm for the shortest path in the graph is original). The proposed numerical methods are mutually compared. A representative sample of adequate problems from pedestrian evacuation is solved. The novelty of the paper is a comparison of the results of different numerical methods for solving the given problem. The article thus deviates from the usual scheme One problem-One method-One result. We have One problem-Four methods-Four results relatively in a good agreement. The choice of method is then up to the user: simplicity versus higher order versus computational time.
The paper is organized as follows. In Section 2 the macroscopic model of pedestrian flow describing the room or open space pedestrian evacuation is formulated. The splitting technique and its numerical realization is explained in Section 3. The physically relevant computations of various evacuation scenarios are presented in Section 4. A conclusion summarizing the results achieved and emphasizing the novelty of the paper is included in Section 5.

Governing Equations
Let us consider the pedestrian flow in the polygonal domain Ω ⊂ R 2 and time interval (0, T) with T > 0. We suppose that the boundary ∂Ω = Γ w ∪ Γ o , where the subscript w states for wall and o for the exit (outflow). The bar in Γ denotes the closure of the set Γ. We further suppose that the subsets Γ w , Γ o ∈ ∂Ω are disjoint and formed by a finite number of connected sets, open with respect to the topology on ∂Ω, where the one-dimensional measure meas(Γ o ) > 0 , i.e., Γ o is not a point. Such a setting is convenient for the modeling of pedestrian evacuation problems in order to predict the pedestrian behavior in panic situations. The system of the PDE's for the vector valued function w is ∂w ∂t Here the unknown variable w = w(x, t) forms a vector w = (w 1 , w 2 , w 3 ) T ≡ (ρ, ρv) T ≡ (ρ, ρv 1 , ρv 2 ) T . We denote by ρ the density of pedestrians in [ped/m 2 ] and by The fluxes f s are defined by The power law for isentropic gases expressing the dependence of the pressure p on the density ρ is used in the definition of the fluxes f s . In our computations the constant p 0 is set p 0 = 1 and the Poisson constant γ = 2. The source term is given by Here V represents the magnitude of the desired velocity in which the pedestrian would like to go, we suppose V = V(ρ) > 0. The notation means that V(x, t) = V(ρ(x, t)) and (x, t) ∈ Ω × (0, T). The definition of V(ρ) can be found in [3]. In our work we use the exponential dependence where α = 7.5, the maximum density ρ max at which is hardly possible to move is set ρ max = 9 and the speed of free flow v max ∈ [1,7]. These parameters setting is based on [4], where experimental investigation of pedestrian behavior is listed. The unit vector µ is the preferred direction of the movement of pedestrians, µ = µ(x, t), (x, t) ∈ Ω × (0, T), and it depends on the distribution of the density ρ in the whole domain Ω, which we express by the functional dependence The function V and the unit vector µ characterize how the speed of pedestrians changes with the density. The parameter τ represents a relaxation time describing how fast pedestrians correct their current velocity v to the desired one Vµ. We set this parametr τ = 0.61 according to the validation of the model with experimental pedestrian movement.

Functional Dependence
To determine the direction µ (desired velocity direction) we have basically two equivalent possibilities how to define functional Ψ in (4).

Integral Formulation
The numerical approximation of φ in (6) is obtained via Dijkstra's algorithm, see original paper [18] or more detailed explanation in [19], and thus we denote it by the subscript D In (6) the minimum is taken over all piecewise regular curves η(σ, x) in Ω, parametrized by the parameter σ ∈ [0, 1] with the property η(0; x) = x, η(1; x) ∈ Γ o . Let us note that the whole distribution of the density at each time instant t has to be given to determine µ.

Differential Formulation
|∇φ( Equation (8) is the so called Eikonal equation of pedestrian flow problem. This equation is completed by the boundary condition (9). The numerical approximation of φ in (8) and (9) is obtained via the Bornemann-Rasch algorithm, see [20], and thus we denote it by the subscript BR.
An alternative approach how to define the functional Ψ was introduced in [21]. The physical meaning of µ as the tangential vector to fastest path to the exit and φ as the shortest time to reach the exit for a given magnitude of velocity V is explained in [8,21].
Let w 0 be a vector function, the initial condition for the system (1)-(4) reads The boundary condition is prescribed where B is appropriate boundary operator. For the boundary condition on the wall we use a normal component of velocity equal to zero v · n = 0.
On the exits we prescribe w 1 = ρ ≈ 0, the space behind the exit is free and empty and pedestrians can leave this space. The more detailed specification of the boundary conditions is described in [8]. Let us note that due to the hyperbolicity of the system the Dirichlet boundary condition cannot be prescribed directly. We solve a linearized initial-boundary value problem instead.

Numerical Solution
The usual way how to solve the hyperbolic systems with the source term is the splitting technique.

Splitting Schemes
We consider two splitting schemes for system (1). Since there are various possibilities of the splitting, a question arises, if the equations are solved "right". This is related to the notion of "verifiation of the model". We answer this question from the numerical point of view. In what follows, we discretize parts of the continuum model by substantially different numerical techniques. Finally we show that the numerical results obtained by the different numerical methods are in a good agreement.
3.1.1. Splitting Scheme I The first splitting scheme for system (1) reads: 3.

Splitting Scheme II
The second splitting scheme for system (1) reads: ∂w ∂t Remember that the first component of w has the physical meaning of density, w 1 = ρ. Above stated splitting Schemes (12)-(14) (Splitting Scheme I) and (15) and (16) (Splitting Scheme II) are employed to numerically solve system (1). Booth schemes are applied in an iterative manner. We employ Finite Volume Method (FV) for the discretization (12)-(14), see Section 3.2, while Finite Volume Method in space and the Discontinuous Galerkin Method in time are used for (15) and (16) (FV-DG), see Section 3.5.
In booth schemes we use triangular mesh T h = {K i } i∈I , where I ⊂ Z + = {0, 1, . . .} is an index set and h > 0. Closed triangles K i , i ∈ I have usual properties from the finite element method. By the discretization a partition 0 = t 0 < t 1 < . . . < t r = T of the time interval [0, T] is meant. We denote by τ m = t m − t m−1 the size of time step between t m−1 and t m and I m = (t m−1 , t m ), m = 1, . . . , r. We denote the set of all vertices of the triangulation T h by V h and the set of all edges of T h by F h .

FV Discretization of the Splitting Scheme I
Since we employ the Finite Volume Method in the Spliting scheme I, the approximate solution is considered as piecewise constant in Ω at the time instant t m Further we denote w m = {w m i } i∈I and derive the FV scheme for w m .

FV Discretization of the Hyperbolic System (12)
Let us assume that w is the exact solution of (12). We integrate (12) over K i , exchange the integration and partial differentiation in the first term and use the Green's theorem with respect to the space coordinates. Then Writing is an index set of neighboring triangles K j to K i (special treatment has to be done if K i has a side on ∂Ω) and approximating where |K i | is the area of K i and H is the so called numerical flux , we derive from (18) the semi-discrete approximation of the hyperbolic system (12) Here |Γ ij | is the length of the side Γ ij . We use the original numerical flux H of Vijayasundaram type proposed in [22]. Its brief description can be found in [8].
For K i having a side on the boundary ∂Ω, the boundary conditions have to be taken into account. On the walls we use inviscid mirror boundary conditions from ( [23], p. 415). Boundary conditions on exit are obtained from solution of a linearized 1D Riemann problem. The formulas for boundary conditions can be found in [8].
If we take Equation (21) for all i ∈ I, we get the vector form of the semi-discrete approximation of (12) d dt where Equation (22) is equipped with the initial condition where w 0 is defined in (10). We denote the approximate solution of (22) at time instant t m by w m , i.e., The following second order Runge-Kutta method can be used for the discretization of Equation (22) where Another possibility is to use the Euler method for the discretization of Equation (22) The initial condition w 0 is set according to (25) The time step τ m is given by the CFL stability condition stemming from the finite volume approach. For the Runge-Kutta method, the time step is halved. Relations (27)-(29) represent the discrete analog to the hyperbolic system (12), respectively. (13) The numerical approximation of µ will be sought in the space of piecewise constant functions on K i at t m , i.e.,

FV Approximation of µ in
Two formulations of the functional Ψ are considered in (4), namely Ψ D and Ψ BR defined in (5)- (9). The subscripts mean that the Dijkstra or Bornemann-Rasch algorithm are used for the computation of Ψ.

FV Approximation of Ψ D
To compute integral in (6) we approximate φ(P) at vertex P ∈ V h by φ P . The approximation φ P is obtained by the minimization process over edges of the triangulation T h only. To this end we approximate the quantity V(·, t m ) on each edge E ∈ F h by a arithmetic mean of V m P s , V m P e at vertices P s , P e of the edge E, see (32). The approximations V m P of point values of V(P, t m ) are obtained from the piecewise constant approximations w m i of w(·, t m )| K i on triangles K i as weighted averages of values on the set of triangles K P having the point P in common. The definition of V and relation (3) is taking into account. Let us note that the first component of w m i is the approximation of the density. Assuming above mentioned leads the integral in (6) can be determined exactly. The meaning of this integral is the time it takes to go from the vertex P to the exit Γ o . The pedestrian path is formed by edges E of the graph given by the triangulation T h and we get Here |E| represents the euclidean norm (length) of the edge E, E is the shortest path in the graph starting from the vertex P and ending in some vertex of Γ o . Since we can define norm · of an edge E the term shortest path means with respect to this norm. According to the physical meaning of quantities in the definiton (33) the shortest means the fastest. The problem of finding shortest path is solved by modified Dijkstra's algorithm where we use norm (33) instead of the length of the edge E. We utilize the fact the graph is not oriented and choose the starting vertex as the vertex of the boundary Γ o and set the ending vertex P ∈ V h . In [24] the multiple starting point Dijkstra's algorithm has been proposed. This algorithm computes the approximations φ P for all vertices P ∈ V h in one run taking an arbitrary vertex lying on Γ o as the starting point. The more detailed description can be found in [8].
From point values φ P , P ∈ V h , the piecewise linear approximation φ h of φ in the space and used for the approximation of µ in (31) For the simplicity of notation the argument in Ψ D,i is written as w m although µ m i depends only on the piecewise constant approximation of the density. Relation (34) expresses the fact that µ m i depends on the approximation of the density on all K j , j ∈ I.

FV Approximation of Ψ BR
We implemented the Bornemann and Rasch algorithm [20] for the solution of the Eikonal Equations (8) and (9). The function φ (called potential) is approximated by the piecewise linear function φ h ∈ X h . The key point of this algorithm is the use of the Hopf-Lax formula in order to obtain the piecewise linear finite element solution of the certain local Dirichlet problem. To this end we construct the approximations V m P of point values of V(P, t m ) from the piecewise constant approximations w m i of w(·, t m )| K i on triangles K i as weighted averages of values on the set of triangles K P having the point P in common.
The sought potential φ h is the fix point of the operator Λ h , where Λ h : X h → X h is the so called Hopf-Lax update function: Here ∂K P is the boundary of the set of volumes having the vertex P in common. The details about the minimization process in Hopf-Lax update Formula (35) and the convergence of the fix point iterations can be found in [20]. Having computed the potential φ h ∈ X h at each vertex of the mesh, we assess the gradient ∇(φ h | K i ) at the centre of gravity of the triangle K i and the quantity Relations (34) and (37) represent the discrete analog of the determination of the desired direction µ in (13). (14) Integrating (14) over K i , we proceed analogously as in the case of the discretization of the hyperbolic system (12). We define the semi-discrete approximation of (14) as

Discretization of System
The notation (19), (23) and (26) and the analogous notation for µ and the source term s is used. The quantity µ m i is obtained from w m using the relations (34) or (37). Equation (38) can be solved numerically by the same Runge-Kutta method as Equation (22). where Another possibility is to use the Euler method for the discretization of Equation (38) Relations (44)-(46) represent the discrete analog of system (14), respectively. The finite volume discretization of the Splitting Scheme I (12)- (14) consists of the passage through discretizations (27)  (3) w := w hyp + τ m S w hyp ; In (1) the following notation is used In (3) the following notation is used The quantity µ hyp i is obtained from w hyp using relations (34) or (37). Analogouslyμ i is obtained from w.

FV-DG Discretization of the Splitting Scheme II
The FV-DG space-time discontinuous Galerkin method was derived in [8]. In the following paragraphs we briefly summarize it. Let q be an integer and P l (0, 1) the space of polynomials of degree l on (0, 1) for l = 0, . . . , q. The approximate solution w h is sought in the finite element space The approximate solution w h ∈ S q h,τ can be written in the form Provided ξ is a function defined in r m=1 I m , we use the notation if the one-sided limits lim where Πw 0 is the piecewise constant projection of the initial condition w 0 on the mesh T h . (15) The desired direction µ in (15) depends on the distribution of ρ(·, t). Two formulations of the functional Ψ are considered in (15), namely Ψ D and Ψ BR defined in (5)- (9). The meaning of subscripts is the same as above.

FV-DG Approximation of µ in
We Having point values V m−1 P the situation is similar to the Section 3.2.3 and with the same notation the approximation φ P of φ(P) is defined To find the shortest path in the graph given by the triangulation T h the modified Dijkstra's algorithm is employed. From the values φ P , P ∈ V h , the piecewise linear function φ h is constructed. This recovery is used to compute gradients of φ h on triangles K ∈ T h . The point approximations of Ψ by weighted average of gradients of φ h is determined where K P is the set of triangles from T h having the vertex P in common and weights in (52) are areas |K| of each triangle K ∈ T h .

FV-DG Approximation of Ψ BR
The point values The passage from w m−1 h | − m−1 to the definition (54) can be written as where the functionals Ψ D,h and Ψ BR,h represent the process described above, see

FV-DG Discretization of (16)
The discretization of (16) is described in detail in [8]. Here we mention the main ideas. Let w be the exact solution of (1). Multiplying (1) by ξ ∈ S q h,τ and integrating over K i × I m , using the Green's theorem with respect to the space coordinates and taking sum over all elements K i ∈ T h we obtain (some manipulation needed and using the fact ξ is constant with respect to the space) where by w we denoted ∂w/∂t. Firstly we discretize the time derivative. Detailed procedure for the treatment of the time derivative can be found in [8]. This procedure is based on the integration by parts. The result of this is the following identity where we use the notation introduced in (49) and ξ = ∂ξ/∂t. Secondly, the integral over ∂K i in (56) is approximated. We use the numerical flux where w| int denote the values of w on ∂K i considered from the interior and exterior of K i , respectively.
The same numerical flux as in (12) is used including the boundary conditions specification.

Resulting FV-DG Scheme
Using both key ingredients, time derivative and numerical flux H together with treatment of the boundary condition for each time interval I m , m = 1, . . . , r, we define the following form where ∂K i \ ∂Ω denotes the part of the boundary of K i ∈ T h lying in the interior of Ω, ∂K i ∩ ∂Ω denotes the part of the boundary of K i ∈ T h lying on ∂Ω and w BC is given by the boundary conditions.

Definition 2.
We say that w h ∈ S q h,τ is the approximate solution of (1) if where µ h ∈ S

Solution of FV-DG Scheme
As the result of the used numerical flux of Vijayasundaram type, we can linearize the form A m h given by (59) as where the form A m L,h is linear with respect to its third argument. The formula for A m L,h contains the linearized numerical flux H and together with a m h (w, ξ) can be found in [8]. Taking into account (48), we write system (60) as which is system of algebraic equations (for each m = 1, . . . , r). Using linearization (61) we approximate w m h by the following sequence of approximations w m,l h , l = 0, 1, 2 . . .
where d m,l h solves This scheme represents a Newton-like method. Here the matrix given by the linearized form A m L,h in (64) is used instead of the Jacobi matrix in the original Newton. For the solution of the linear system (64) the GMRES method with block ILU(0) preconditioner is used.
To improve the convergence the adaptive damping factor λ ∈ (0, 1] in (63) is used w m,l+1 h = w m,l h + λd m,l h . More details about adaptive dumping factor can be found in [25,26]. In [8] we discuss stopping criteria for the used iterative solver.

Results
We present the comparative study of the proposed Splitting schemes I and II from the numerical point of view. The described methods are compared on domains with different geometries. Several meshes grading from coarse to fine are used. Let us stress that we have tested two different numerical methods on the representative sample of computational domains. Moreover, the software for the Finite Volume Method was written by P. Kubera (the first author), while the software for the Discontinuous Galerkin Method was written by V. Dolejší and it was used as the black box with his permission. The detailed description of Dolejší's method can be found in [8] and generally in [23]. In each splitting scheme I and II the two different approaches for the determination of µ, namely Bornemann-Rasch and Dijkstra's algorithms, were used. So we gained a more general insight on the behavior of different methods for pedestrian movement. See also the further comparison with the new Cellular Automaton model proposed recently by the authors in [15].
Quantities we are interested in are the distribution of the density ρ (the dimension being pedestrian per square meter in the distinction of the fluid dynamics density in kg/m 3 ) over the domain and the evacuation time in seconds T evac , i.e., the time when the domain is considered to be empty of any people. See the first paragraph in Section 3 for the details how the results have been computed. For each of test cases we use two methods for the determination of µ in (4). For the numerical solution of φ in the integral Formulation (6) the abbreviation DA is used and it suggests the Dijkstra's algorithm. For the numerical solution of φ in the differential Formulation (8) the abbreviation BR indicates the Bornemann-Rasch algorithm.

Let us define the quantity M(t m )
giving the number of pedestrians in the domain Ω at the time instant t m . The quantity M(t m ) is approximated numerically using the computed piecewise constant approximation of the density at time t m . We consider the domain Ω to be empty if the quantity M(t m ) ≤ #PEDS, where #PEDS is a prescribed number. We set #PEDS = 2 for all methods.

Domain with One Obstacle
The first example is the evacuation of the rectangular domain 40 × 10 m with one circular obstacle, see Figure 1  The results concerning the evacuation time T evac are summarized in Table 1. We see that the difference in T evac is less than two seconds.  In Figure 2, the evolution in time t = 10, 15, 20 s of the pedestrian density is drawn. The intuitive movement of pedestrians towards the right hand side can be observed. In each row of figures, the FV and the FV-DG methods are compared. We see all methods are in a good agreement. In Figure 3 the time evolution of the total number of pedestrians M(t m ) is shown.

Domain with Three Obstacles
Here we present results for the simulation of the evacuation of the rectangular domain 40 × 10 m with three circular obstacles, see Figure 1 (bottom). We suppose boundary and initial conditions and also parameters of the model to be the same as for one circular obstacle above.
Computational results are shown in Table 2. We can see the difference in T evac being smaller than one second for all methods on all tested meshes. The evolution of the pedestrian density in time is depicted in Figure 4. The intuitive behaviour with the high-density domain in front of the obstacles can be seen for all methods. The dependence of the total number of pedestrians M(t m ) for all tested methods is shown in Figure 5. The curves almost coincide (left). The small difference can be seen on the close view (right). We see that all methods are in a good agreement.

H-Shape Domain Evacuation
Further we present an evacuation through a narrowing followed by an empty space with the non-symmetric exit, see Figure 6. The so-called H-shape domain represents an example of the non-convex domain. The convexity of the domain does not play any role in the design of the proposed numerical methods. The geometry is defined The outflow boundary condition is prescribed on and the wall boundary condition on Γ w = ∂Ω \ Γ o . We prescribe the initial condition 25], 0 elsewhere, and v 0 = 0.  Table 3 shows the computational results obtained by the FV and by the FV-DG methods. See the description of table, the use of the DA or BR algorithm in the pedestrian flow model (both FV and FV-DG) causes more significant differences in the evacuation time T evac . On the other hand the FV and FV-DG model with the same type of the computation of the desired direction of motion (i.e., DA or BR) produce the T BR evac and T DA evac being in a good agreement. The dependence of the pedestrian evacuation model on the use of the DA or BR algorithms has to be further studied. In Figure 7 the density evolution for selected time instants is depicted. Here the local increase of the pedestrian density in the left part of the narrowing and also close to the exit of the domain is noticeable. In Figure 8 the evolution of number of pedestrian in the domain is shown. The main differences are due to the use of the DA or BR algorithms.

T-Shape Domain Evacuation
Here we present numerical example for the multiple exits non-convex domain, the socalled T-shape domain, see Figure 9 for details. The  Table 4 shows the computational results for the T-shape domain. There is a small difference in T evac between FV and FV-DG but the results are still in a good agreement from the practical point of view. The evolution of number of pedestrians in T-shape domain is shown in Figure 10. The pedestrian density distribution over the domain for different time instants is depicted in Figure 11.   The Figure 12 shows the trajectories of "selected" pedestrians. Trajectories of the pedestrians are given by the solution of the differential equation with initial condition x(0) = x 0 . The initial condition arranges the pedestrian selection. The Equation (69) is solved by the forward Euler method. Here we use the piecewise constant approximation of v(x, t) obtained from the numerical solution w m h and w h (see Definitions 1 and 2). The trajectories in Figure 12 show us the influence of the density distribution on the movement of pedestrians and their trajectory, since µ depends on the density.

Conclusions
The paper deals with the comparative study of two splitting methods-Splitting Method I and Splitting Method II. They were described and numerically tested for solving several pedestrian flow problems. The FV approach was compared with FV-DG approach. The former use the piece-wise constant approximation together with the explicit time discretization and is advantageous from the implementation point of view. The latter is based on the discontinuous piece-wise polynomial approximation and uses the implicit discretization by space time discontinuous Galerkin method. The FV-DG approach requires the solution of the system of nonlinear algebraic equations, but on the other hand the implicit time discretization together with the high polynomial degree enable us to use the longer time steps. Further the FV-DG approach allow us to use different space grids on different time levels which is useful for adaptive methods.
The alternative method to that proposed by Bornemann and Rasch for the computation of the desired direction of move µ in (4) based on the Dijkstra's algorithm was proposed and tested. This method can be considered as a competitive method to the Bornemann-Rasch algorithm for the solution of the Eikonal equation. The key difference is in leaving out the differential formulation given by the Eikonal equation (8) and replacing its iterative solution by the modification of the Dijkstra's algorithm for the shortest path in a graph. The curvilinear integral in (6) is minimized in the class of piece-wise linear curves given by edges of the triangulation. Compared with the Bornemann-Rasch algorithm (BR) the proposed method has several different properties. Firstly, the method is tolerance free in the distinction from the iterative BR algorithm, where the stopping criterion must be prescribed. Secondly, in the proposed method the computation of the gradient of φ in vertex P can be avoided by using the direction of the edge which lies on the shortest path from P to the outflow.
The special advantage of presented simulations is that the system of partial differential equations is solved by two different software systems with the results being in the good agreement. This is related to the key word "verification" of the model, i.e., to the question, if the equations are solved "right". Another step is the question, if we solve the "right" equations, which is related to the notion "validation" of the model. This suppose the comparison with the experimental data. Two ingredients of the proposed methods are original. We use the numerical flux which is the original modification of the Vijayasundaram numerical flux for the Euler equations. Despite the fact that in the distinction from the Euler fluxes the fluxes for the pedestrian flow equations are not homogeneous, we were able to propose the numerical flux of the Vijayasundaram type which is easy to be coded. For details see [8,22]. Further we are able to avoid the numerical solution of the Eikonal equation. This is related to the notion of the fastest path in a graph. The details can be found in [8,27]. We use the triangular finite volumes and the numerical flux mentioned above. In [3] the dual finite volumes are used and the HLL numerical is flux applied. The results can be mutually compared.