Optimal number of Schur subdomains: Application to semi-implicit finite volume discretization of semilinear reaction diffusion problem

The purpose of this paper is to establish a new numerical approach to solve, in two dimensions, a semilinear reaction diffusion equation combining finite volume method and Schur complement method. We applied our method for q = 2 non-overlapping subdomains and then we generalized in the case of several subdomains (q≥2). A large number of numerical test cases shows the efficiency and the good accuracy of the proposed approach in terms of the CPU time and the order of the error, when increasing the number of subdomains, without using the parallel computing. After several variations of the number of subdomains and the mesh grid, we remark two significant results. On the one hand, the increase related to the number of subdomains does not affect the order of the error, on the other hand, for each mesh grid when we augment the number of subdomains, the CPU time reaches the minimum for a specific number of subdomains. In order to have the minimum CPU time, we resorted to a statistical study between the optimal number of subdomains and the mesh grid.

Abstract. The purpose of this paper is to establish a new numerical approach to solve, in two dimensions, a semilinear reaction diffusion equation combining finite volume method and Schur complement method. We applied our method for q = 2 non-overlapping subdomains and then we generalized in the case of several subdomains (q ≥ 2). A large number of numerical test cases shows the efficiency and the good accuracy of the proposed approach in terms of the CPU time and the order of the error, when increasing the number of subdomains, without using the parallel computing. After several variations of the number of subdomains and the mesh grid, we remark two significant results. On the one hand, the increase related to the number of subdomains does not affect the order of the error, on the other hand, for each mesh grid when we augment the number of subdomains, the CPU time reaches the minimum for a specific number of subdomains. In order to have the minimum CPU time, we resorted to a statistical study between the optimal number of subdomains and the mesh grid.
applications in physics in [3,16]. A general semilinear reaction diffusion equation can be expressed as ∂ t c(t, x) − ∇(D(x)∇c(t, x)) + f (c) = g(t, x); t > 0, x ∈ R N , N ∈ N * , where D is the diffusion coefficient, g(t, x) is a given source term and f (c) is the nonlinear reaction term of the equation.
Several works have treated the resolution of the nonlinear reaction diffusion equations, we refer for instance to [4] where a two-grid finite volume element methods for semilinear parabolic problems has been presented and studied by Chen and Liu and [7] where the Fernandes and Fairweather described an alternating direction implicit orthogonal spline collocation method for nonlinear reaction diffusion systems. The discretization of these problems leads to large systems, hence the use of the domain decomposition methods (DDM) is needed. The first models of these methods have been established by Schwarz [14,5], and after that, several variants of the Schwarz method have appeared [9,6,13], for example, in [1] the authors used the non-overlapping domain decomposition algorithms of Schwarz waveform relaxation type to resolve the semilinear reaction-diffusion equation.
In the first part of this work, we establish a semi-implicit finite volume scheme (FV). The linear part is discretised by the Euler implicit scheme and the non linear reaction term is integrated explicitly. To improve the calculation cost, we propose a new approach (FV-SC) coupling the last scheme and the algebraic Schur complement technique (see, for instance, [2]), this method was applied (developed and studied) for a linear parabolic advection diffusion problem (see, for example, [10]). Several numerical experiments show the interest of the proposed algorithm in terms of accuracy, robustness and efficiency. All results have been validated using a numerical test cases and also by using the free finite elements software Freefem++. In the second part, we discuss the obtained numerical results and we show experimentally the accuracy and the efficiency of the proposed method, measured both in order of error, and in CPU time. Based on the obtained results, and using a statistical study, our challenge is to find the relation between the optimal number of subdomains and the mesh grid, in order to have the minimum CPU time. Let's emphasize that up to our knowledge, no work dealing with this study have been presented in the literature.
2. Problem description. We are interested with the following semilinear reaction diffusion equation: where Ω ⊆ R 2 , T > 0 and the diffusion coefficient ν > 0. The function f is a nonlinearity which belongs to C 2 (R), and satisfies f (0) = 0. The initial data c 0 is supposed to be defined in H 2 (R 2 ). In the case where Ω is a bounded subset of R 2 , we need to specify a boundary condition on ∂Ω. It can be, for example, a Dirichlet condition c(x, t) = c D (x, t) on ∂Ω × (0, T ). It is easy to verify that the problem (1) has a weak solution c which belongs to L 2 (0, T ; H 1 (R 2 )) ∩ C([0, T ]; L 2 (R 2 )), such that f (c) ∈ L 2 (0, T ; L 2 (R 2 )), satisfying where ( , ) denotes the inner product in L 2 (R 2 ). Let us recall the following result concerning the well-posedness of the Cauchy problem (1).
Proof. It is an easy task. One can see [3].
3. Finite volume approach. The finite volume approach consists in dividing the domain of calculation Ω (bounded subset of R 2 ) into a finite number of control volumes (CVs) We are interested in the so called cell centered FV approach, where the discrete unknowns are located at the center of CVs. For a general CV, we use the notation of the distinguished points (midpoint, midpoints of faces) and the unit normal vectors according to the notation as indicated in figure 1 (right). We denote the midpoints of neighboring CVs with capital letters W , S, etc (see figure 1 left), we adopt the notations given in [10].
where S a (a = e, n, w, s) are the four faces of volume V P (figure 1), n a are the unit normal vectors to the face S a . The linear part of the equation is discretized using implicit finite volume scheme and the nonlinear term f (c) is integrated in the scheme explicitly. We obtain where µ(V P ) is the volume of cell V P , F n+1 P = Sa (−ν∇c n+1 )n a dS a denotes the diffusion fluxes through the CV faces, and c 0 • For discretization of diffusion term, we have considered a centred difference scheme, for example,  (3) can be arranged as follows: ψ P is a constant depending on the boundary and initial conditions. The coefficients a I (I=P , E, W , N or S) are defined as follows: Finally, the numerical scheme is expressed as the following linear system: C n , f (C n ) and Ψ are the vectors of c n P , f (c n P ) and ψ P (respectively).

Non-overlapping domain decomposition.
To reduce storage and calculation time, the domain of calculation is usually decomposed into numerous subdomains. The domain Ω is decomposed into multidomain non-overlapping strip figure 2). Let's denote by Γ ij the interface between Ω i and Ω j , with Γ = ∪ ij Γ ij , and by n i the normal direction (oriented outward) on Γ ij for i = 1,· · · , q − 1 and j = i + 1. For simplicity of notation we set n = n i .

Figure 2. Non-overlapping strip decomposition
Considering a rectangular mesh of Ω, each subdomain Ω i is partitioned into n i (i = 1,· · · ,q) cells in X direction and m cells in Y direction (figure 3).
The problem (1) can then be expressed as: The last two interface conditions are known as transmission condition on Γ ij . The decomposed problem (5) is discretized on each subdomain Ω i , i = 1, · · · , q using the finite volume scheme described in Section 3.
For the interface conditions, the centered differences scheme has been used. The following system, for i = 1,· · · , q − 1 and j = i + 1, is obtained: where σ i = e i and σ j = w j if V P i ∩ Γ ij = ∅ (i = 1, · · · , q − 1), σ i = E i and σ j = W j otherwise, and ψ P i are constants depending on initial and boundary conditions on Ω i (i = 1,· · · ,q).
5. Schur complement. The Schur complement method is a kind of non-overlapping DDM, that obtains solutions using interfacial system constructed by static condensation. This interfacial system, known as the Schur complement system, is easily formulated as shown below. The methods based on Schur complement exists in two versions. The first one uses the Steklov Poincaré operator and the second one is an algebraic version.
For example in [10] one finds presentation of these methods used in the context of a finite volume method (for linear advection-diffusion equation) and in [14,15,19], a presentation of these methods used in the context of a finite element method.
In this work, an algebraic version of Schur complement method has been used. Let C n+1 i and C n+1 Γ denote the vector of the unknowns of Ω i (i = 1,· · · ,q) and Γ at time t n+1 (respectively), so the decomposed problem (6) can be written in the 26 HASSAN BELHADJ, MOHAMED FIHRI, SAMIR KHALLOUQ AND NABILA NAGID following matrix form: with A i , A iΓ , A Γi and A ΓΓ (i = 1,· · · ,q) describe respectively (a), (b), (c) and (d) of system (6). f (C n i ) and Ψ i are respectively the vectors of f (c n P i ) and ψ P i (i = 1, · · · , q −1). A i presents the coupling of the unknowns in Ω i , A ΓΓ is related to the unknowns on the interface, A Γi and A iΓ represent the coupling of the unknowns of each subdomain Ω i with those of the interface Γ ii+1 , for i = 1, · · · , q − 1.
The system (7) can be solved formally by block Gaussian elimination. Eliminating C n+1 i (i = 1,· · · ,q) in the system (7) yields to the following reduced linear system for C n+1 The matrix S is the Schur complement matrix. After calculating C n+1 Γ , C n+1 i can be obtained immediately and independently (in parallel) by solving: Γ , (i = 1, · · · , q). In short: • FV-SC algorithm: We assume that Ω is decomposed into q ≥ 2 subdomains Ω i=1,...,q . For each time iteration n, the following steps are performed: 1. Construct blocks A i , A iΓ and A Γi , 2. Calculate S and χ Γ by removing all internal unknowns C n+1 i=1,...,q 6. Numerical simulations. In this section, the spatial domain is the square Ω = (−1, 1)×(0, 2), decomposed into q ≤ 9 rectangular subdomains. The time interval is (0, 1) and the diffusion coefficient ν equals 1. We test here the case of the nonlinear function f (c) = c 3 (studied e.g. in [1]) with the Dirichlet boundary conditions (e.g. c D = 2) and a given constant initial data (e.g. c 0 = 2). All simulations are performed on a structured and conforming square mesh of the space domain (see figure 1 and figure 3). Different values of the spatial mesh step h are considered, and the time step ∆t is chosen such that ∆t = h. Note that for this problem, there is no analytical known solution. The numerical experiments are compared to the finite element method using the free software Freefem++ 3.20 dedicated to this kind of methods.  Several test-cases have been performed to validate the FV mono-domain (FV-1D) algorithm and FV coupled to the Schur Complement method (FV-SC). The coupled FV-SC algorithm was tested for q subdomain decomposition (resp. q = 2, . . . , 9). For all the test-cases, we have considered different values of h and of time steps t. Let use the following notations: • FE-Pi is the finite element method with the basis functions in the set of polynomials of degree i (=1,2), • h F E is the spatial mesh step for FE method, • h F V is the spatial mesh step for FV method, • the spatial mesh step is denoted by h when h F E = h F V . We can see in figure 4 (respectively figure 5) the approximation solutions of FE (left) and FV-1D (right) in 3 dimensions (respectively 2 dimensions). Figure 6, shows the point-wise differences between FE-Pi and FV-1D, as for figure 7, it shows the point-wise differences between FE-Pi and FV-SC. In the latter figures we observe that the approximation solutions by FV-1D and FE-Pi have the same shape, furthermore the order of the point-wise differences FE(Pi)-FV(1D) and FE(Pi)-FV(SC) are 10 −4 . In figure 8, we have plotted the L ∞ (Ω) error history between FE-Pi (i=1,2) and FV-1D (left) (FV-SC (right), respectively), for h = 0.0138 and t = 1. We remark that, for both situations, L ∞ -error is small when the functions basis are in P2 which ensures the good accuracy of our code.   To verify the efficiency of the proposed algorithm, we also have computed the CPU time of calculation, when increasing the subdomain number. We have presented the corresponding results in table 2 and figure 11. Note that the numerical experiments was performed with different values of the mesh step h F V (h F V = 0.1666, 0.0556, 0.0333, 0.0062) and with different number of subdomains (q ≤ 9). From these numerical results, one can see an important advantage concerning the computation time of the FV-SC algorithm. Indeed, we have obtained a gain of almost: • 50% CPU time for h = 0.0138 and 4 subdomains, • 45% CPU time for h = 0.0069 and 8 subdomains.  7. Study of optimal number of subdomains. According to the numerical results of Section 6, it was noted that by increasing the number of subdomains, the CPU time reaches the minimum, for a specific number of subdomains. We have presented the corresponding results in table 3. Based on the obtained results, and using a statistical study, we aim to find the relation between the optimal number of subdomains (SD) and the mesh grid (h), in order to have the minimum CPU time.
A linear fit (linear regression), or intrinsically linear (logarithmic, power, exponential, etc.) is studied. In order to choose the best model adjusting the observations, the correlation coefficient and the coefficient of determination are calculated [8,18]. The correlation coefficient (denoted by R), developed by Karl Pearson, measures the strength and the direction of a linear relationship between the two variables, we say that we have a strong positive (respectively, negative) linear correlation if R is close to +1 (respectively, −1). The coefficient of determination is equal to the square of the correlation coefficient (R 2 ). It represents the percent of the data that is the closest to the line of best fit. It is a measure of how well the regression line represents the data. Based on the observed data, we found that the relationship between the two variables is well described using the power model, and takes the form SD = 0, 4645343263 × h −0,5235932037 .
Graphically, this function is illustrated in figure 12. By applying the logarithm, In this case, we found that the correlation coefficient is equal to R = −0.99, the negative value indicates that as values for ln h decrease, values for ln SD increase. The coefficient of determination is equal to R 2 = 0.98, this means that 98% of the total variation in ln SD can be explained by the linear relationship between ln h and ln SD, this is a quasi-perfect representation. The 2% of the total variation in ln SD remains unexplained.

Conclusion.
A new approach coupling semi-implicit FV and Schur complement methods applied to a semilinear reaction diffusion equation, on 2D structured and matching mesh, is presented. The algorithm applied to non-overlapping multidomains decomposition has the proprieties of stability and convergence. On the other hand, we preserve a good accuracy of calculation with the increasing of the subdomains number comparing with global calculation FV-1Dom. Furthermore, we have a better CPU time when increasing the subdomains number, for example: • gain of almost 50% CPU time for h = 0.0138 and 4 subdomains, • gain of almost 45% CPU time for h = 0.0069 and 8 subdomains. A statistical study allows us to find a relation between the optimal number of subdomains and the mesh grid, in order to have the minimum CPU time.
As perspectives, we aim to: • use unstructured and matching mesh for complex applications, • use unstructured and nonmatching mesh for complex applications, • make a theoretical study of the relationship between the optimal number of subdomains and the mesh grid.