NONLINEAR ELLIPTIC PROBLEMS WITH THE METHOD OF FINITE VOLUMES

We present a finite volume discretization of the nonlinear elliptic problems. The discretization results in a nonlinear algebraic system of equations. A Newton-Krylov algorithm is also presented for solving the system of nonlinear algebraic equations. Numerically solving nonlinear partial differential equations consists of discretizing the nonlinear partial differential equation and then solving the formed nonlinear system of equations. We demonstrate the convergence of the discretization scheme and also the convergence of the Newton solver through a variety of practical numerical examples.


Introduction
Nonlinear elliptic equations arise in many applications in many fields, so solving such systems is important.Solving nonlinear partial differential equations consists of discretizing the partial differential equations and solving the formed nonlinear algebraic system of equations.Work has been done on numerically solving nonlinear elliptic partial differential equations (PDEs).For example, Schwarz alternating methods (see [14] and references therein), multigrid methods [5], and preconditioned FFT [12].In this work, we explore the convergence of the discretization method, and also the convergence behaviour of the Newton-Krylov method for solving the nonlinear algebraic equations [7,8].Let us consider the following nonlinear elliptic problem: p(x, y) = p D on ∂Ω D , g(x, y) = −K∇p on ∂Ω N . (1.2) The above problem captures the fundamental features of the Poisson-Boltzmann equation arising in molecular biophysics (see [2][3][4][5][6]).Here, Ω is a domain in R 2 , the source function s is assumed to be in L 2 (Ω), and the diagonal tensor coefficient K(x, y) is positive definite and piecewise constant.K is allowed to be discontinuous in space.In biophysics literature, the medium properties K is referred to as the permittivity [2][3][4][5].It takes the values of the appropriate dielectric constants in the different regions of the domain Ω.In (1.1) and (1.2), ∂Ω D and ∂Ω N represent the Dirichlet and Neumann parts of the boundary, respectively.f (p) represents the nonlinear part of the problem.Equations (1.1) and (1.2) model a wide variety of processes with practical applications, some examples are pattern formation in biology, viscous fluid flow phenomena, chemical reactions, biomolecule electrostatics, and crystal growth.This paper presents the two-point finite volume discretization (2P-FVM) [9,10] of the nonlinear problem (1.1) and (1.2) on the rectangular meshes.An implementation of boundary conditions is also presented.Several numerical examples are reported for showing the convergence of the finite volume discretization scheme.A Newton-Krylov algorithm is also mentioned for solving the system of nonlinear equations formed by the discretization scheme.Convergence of the Newton method is demonstrated through numerical work.For higher-order finite volume discretization of linear problems, the interested reader are referred to [1,11].Handling complex geometries is a difficult task; radial basis functions (RBFs) are a new numerical method that can offer very high accuracy even on complicated domains [13].RBFs can be promising in solving nonlinear partial differential equations.
An outline of the paper is as follows.In Section 2, finite volume discretization of the nonlinear elliptic equations is presented.An implementation of Neumann and Dirichlet boundary conditions is also reported.Section 3 presents a Newton-Krylov algorithm for solving the system of nonlinear equations.Numerical work is reported in Section 4. Finally Section 5 concludes the paper.

Two-point finite volume discretization of the nonlinear elliptic problem
For solving partial differential equations (PDEs) on a domain by numerical methods such as the 2P-FVM, the domain is divided into smaller good quality elements (meshing of the domain).These elements are called finite volumes or cells.The degrees of freedom (DOFs) for the 2P-FVM lie at the cell centers.Thus, each finite volume in the mesh gives rise to an algebraic nonlinear equation corresponding to (1.1).A residual form of (1.1) is − div(Kgrad p) + f (p) − s(x, y) = 0. Integrating it over one of the finite volumes V in the mesh and using the Gauss divergence theorem lead to where n is the outward unit normal on the boundary (∂V ) of the cell V .Let us assume that finite volumes V are rectangular or quadrilateral in shape.Thus, the boundary of these finite volumes consists of four segments and the above equation can be written as Sanjay Kumar Khattri 3 The term ∂Vi −K∇p • n is referred to as the flux through the edge ∂V i .Let us denote it by Ᏺ i .Thus, (2.2) can be written as (2. 3) The degrees of freedom for the 2P-FVM lie at the cell centers, so the scalar variable p is assumed constant in each cell.Thus, the integral V f (p) is approximated as V f (p) ≈ f (p)V .Here, V is the area of the cell.Thus, the above equation can be written as For evaluating the integral V s, we use the Newton-Cotes formulas.Each finite volume in the mesh results in the nonlinear equation (2.4). Figure 2.1 shows the equation for a cell.Collecting all such nonlinear equations will result in a discrete system of nonlinear equations A(p h ) = 0. Now let us consider computing Ᏺ i in (2.4).
where Φ 12 is a scalar and is given as Here, l is the length of the common edge AB. h 1 is the perpendicular distance from the center of the cell P on the edge AB.Similarly h 2 is defined.Let there be total n cells in the mesh (DOF = n).Each cell in the mesh provides a nonlinear equation (2.4).Collecting all such equations results in a system of nonlinear equations given as The above nonlinear system of equations can be solved by a Newton-Krylov method [8].
In Section 3, a Newton-Krylov algorithm is presented for solving the nonlinear system A(p) = 0.

Implementation of boundary conditions.
In the case of finite volume discretization, every finite volume in the mesh results in a nonlinear discrete equation ( 2 Here, Ω is the area of the triangle, and n i is the normal vector on the edge opposite to the vertex i.The magnitude of the vector n i is equal to the length of the edge.Let the unknown potential at the center of the boundary cell 3 be p 3 .The potential gradient inside the boundary triangle can be approximated by the expression (2.8).Thus, the flux through the boundary edge 12 is F 12 = −(K∇p) • n 3 .Let the outward normal vector on the edge (see Figure 2.3(b)) opposite to the vertex i be n i = (nx i ,ny i ) t .The vector n i is pointing away from the node i and the magnitude of the vector is equal to the length of the edge.Let the property K of the boundary cell 3 be

Newton-Krylov algorithm
The nonlinear algebraic system of equations A(p) can be expanded by the Taylor's series around an initial guess p 0 as where J is the Jacobian matrix, HOT exists for higher-order terms, J(p 0 ) is the value of the Jacobian matrix at the initial guess p 0 , and the difference vector Δp = p − p 0 .The Jacobian J is an n × n (n is the DOF) linear system.The Jacobian J is given as Here, p 1 , p 2 ,..., p n are the potential associated with the n cells in the mesh.The Jacobian is symmetric, and it will be positive definite and diagonal dominant for positive nonlinearities, that is, for f > 0. Setting (3.1) equal to zero and neglecting the higher-order terms results in a basis for the Newton algorithm The above linear system is a basis for the Newton algorithm for finding the zeros of the nonlinear function A(p).The linear system (3.3) is solved by the conjugate gradient solver [15, chapter 5].The Newton iteration for solving A(p) = 0 is A Newton-Krylov iteration for solving A(p) = 0 is given by Algorithm 3.1.In Algorithm 3.1, • L2 denotes the discrete L 2 norm and max iter is the maximum allowed Newton's iterations.It is interesting to note the stopping criteria in Algorithm 3.1, we are using three stopping criteria in the algorithm.Apart from the maximum allowed iterations, we are using L 2 norm of residual vector ( F(p) L2 ) and also L 2 norm of difference in scalar potential vector ( Δp L2 ) as stopping criteria for the algorithm.Generally in the literature maximum allowed iterations and the residual vector are used as stopping criteria.If the Jacobian is singular, then the residual vector alone cannot provide a robust stopping criterion.We have implemented the algorithm in the C ++ language.

Numerical examples
Let p be the exact solution vector and let p h be the finite volume solution vector on a mesh.Let us further assume that p k denotes the exact solution at the center of the cell k and p k h denotes the discrete solution by the finite volume approximation for the same location.The error in the L ∞ norm is defined as and error in the L 2 norm is defined as Here, Ω k is the area of the finite volume k in the mesh.
For solving the Jacobian system, we use the conjugate gradient (CG) solver with the ILU preconditioner.For all the numerical examples, our stopping criteria in Algorithm 3.1 are tol = 10 × 10 −20 and max iter = 20.The tolerance for the CG solver is 1 × 10 −10 .During numerical experiments, we observed that the algorithm does not converge for all initial guess of p as it is expected with Newton's iteration.report the convergence behaviour of the 2P-FVM method for the nonlinear elliptic equation.It can be seen in these figures that the method is converging as the mesh is refined.In the L 2 and L ∞ norms, the method is super convergent.
Example 4.2.In this experiment, we solve (4.5) in Ω = (0,1) × (0,1).Let the exact solution be p = (x 2 − x 2 )sin(3π y).Solution inside the domain is enforced by the Dirichlet boundary condition and the source term.We assume zero initial guess for the Newton-Krylov algorithm:       There is a singularity in the solution at origin as it can be expected for the interface problems [10].Figure 4.12 shows the convergence of the Newton-Krylov Algorithm 3.1.Again we see that 5-6 Newton's iterations are sufficient.It should be noted here that problems with discontinuous can produce badly conditioned Jacobian systems: p(x, y) = x(x − 1)y(y − 1) on ∂Ω.(4.8)  .In this example, the source function exhibits a huge variation inside the domain unlike the previous example (source is zero).The source is f = 2.0y(y − 1) + 2.0x(x − 1) − (x(x − 1)y(y − 1))exp(x(x − 1)y(y − 1)).In this experiment, the exact solution is not known.The domain Ω is divided into four equal subdomains (see Figure 4.9) based on the medium properties .Further assume 1 = 1.0I,

Conclusions
We have presented the two-point finite volume discretization of the nonlinear elliptic problems.An implementation of Dirichlet and Neumann boundary conditions is also mentioned.A Newton-Krylov algorithm for solving the system of nonlinear equations is given.Reported numerical work validates the convergence of the discretization scheme for nonlinear elliptic problems.Convergence of the Newton-Krylov method is also reported.

Figure 2 .
Figure 2.1.A nonlinear discrete equation for a finite volume.Here, Ᏺ i with i = 1,...,4 are the fluxes through the cell boundaries.

Figure 2 .
2 shows two cells P and E. Let us compute the flux across the common edge AB of the cells.Let the K of the cells P and E be K P I and K E I. Here, I is the identity matrix.Flux across the edge AB by the 2P-FVM is given as[10]

Figure 2 . 3 .
Figure 2.3.Implementation of the Dirichlet boundary condition.(a) A 3 × 3 mesh.Pressure is specified at the boundary points 1 and 2. Flux (Ᏺ 12 ) through the edge 12 is expressed as a linear combination of the potentials at the locations 1, 2, and 3; see (2.10).(b) Boundary triangle.Here, n i with i = 1,...,3 are the normal vectors on the edges.

Figure 2 .
3(a)  shows a 3 × 3 mesh.Let the pressure be specified at the boundary points 1 and 2. For applying the finite volume formulation (2.4) to the boundary cell 3, we have to compute the flux (Ᏺ 12 ) through the boundary edge 12.For computing the flux, let us form a boundary triangle 123 as shown in Figure 2.3(b).

Example 4 . 1 . 8 Finite
Let us solve the following problem:− p + γpe p = f in Ω, (4.3) p(x, y) = p D on ∂Ω.(4.4)AssumeK = I and the nonlinear function f = γpe p in (1.1).Let the exact solution be p = x(x − 1)y(y − 1) and γ = 1.0.Let the domain of definition be Ω = (0,1) × (0,1).Solution inside the domain is enforced by the source term and the Dirichlet boundary condition.The initial guess for starting the Newton-Krylov algorithm is the zero vector.For understanding the convergence behaviour of the finite volume method, we performed
, 4.2, 4.3, and 4.4 report the outcome of our numerical work.Figures 4.1

Figure 4 . 3 .
Figure 4.3.Example 4.1: convergence of the discretization scheme.L ∞ error versus the degrees of freedom in the mesh.We are observing p − p h L∞ ≈ Ch 2.90 .Here h denotes the size of the smallest edge in the mesh.

Figure 4 . 4 .
Figure 4.4.Example 4.1: convergence of the discretization scheme.L 2 error versus the degrees of freedom in the mesh.We are observing p − p h L2 ≈ Ch 1.99 .Here h denotes the size of the smallest edge in the mesh.

3
and 4.4 report the convergence of the finite volume method.It is interesting to notice in Figures 4.1 and 4.2 that for the three meshes, 5-6 Newton iterations are sufficient for reducing the residuals A(p) L2 and Δp L2 to about 1 × 10 −15 .The convergence of A(p) L2 and Δp L2 is asymptotically quadratic.Figures 4.3

Figure 4 . 7 .
Figure 4.7.Example 4.2: convergence of the discretization scheme.L 2 error versus the degrees of freedom in the mesh.We are observing p − p h L2 ≈ Ch 2.00 .

Figure 4 .
Figure 4.11 is showing the surface plot of the discrete solution on a 64 × 64 mesh.There is a singularity in the solution at origin as it can be expected for the interface problems[10].Figure4.12 shows the convergence of the Newton-Krylov Algorithm 3.1.Again we see that 5-6 Newton's iterations are sufficient.It should be noted here that problems with discontinuous can produce badly conditioned Jacobian systems:
Figure 2.2.Flux Ᏺ AB through the interface AB shared by the cells P and E; see (2.5).
Ω is the absolute value of the area of the boundary triangle 123.6 Finite volume for nonlinear elliptic problems 2.1.2.Flux boundary condition.Implementation of Neumann or flux boundary condition is even simpler.Flux across a boundary edge will be added with the source term V s in the discrete nonlinear equation (2.4).