Conservative Enforcing Positivity Algorithm for Nonlinear Diffusion Equations

1. Conservative Enforcing Positivity Algorithm for Nonlinear Diffusion Equations Guangwei Yuan Institute of Applied Physics and Computational Mathematics, Beijing, China Cooperated with Yanzhong Yao and Fujun Cao

2. Content • Background • Algorithm: Enforcing Positivity with Conservation • Numerical results • Conclusions

3. 1. Background • Nonlinear diffusion equations occur in many different physical systems and applications, such as porous flow, biological systems, inertia confined fusion, and astrophysics. • Diffusion equations possess a remarkable property, i.e., maximum principle, which provides a useful criterion in the study of numerical methods. • The preservation of the maximum principle in numerical solution is one of the key requirements for discrete schemes.

4. If Then A special case is that the solution of diffusion problem is nonnegative if the initial and boundary data and source are nonnegative, which is so called positivity (i.e. nonnegativity) preserving. Our aim is to ensure that the property of positivity-preserving holds for numerical solution.

5. Cell centered finite volumeschemes • 2D rectangle mesh, and scalar diffusion coefficient: Five-point scheme is positivity preserving. • 2D distortedmeshes, or anisotropic diffusion： Many existing discrete schemes are not positivity preserving. Nine point schemes on quadrilateral meshes: Kershaw, Diamond, SOM (MFD), MPFA

6. Many schemes are not positivity preserving:spurious negative solution, oscillation on distorted meshes Diamond MPFA The diffusion coefficient is defined only for nonnegative (e.g., temperature). The calculation stops for !

7. Some known methods • New schemes: positivity preserving • The classical five point scheme is positivity preserving but inaccurate on distorted meshes. • It is difficult to construct a positivity preserving scheme without accuracy loss. Nonlinear monotone schemes • C.LePotier,C.R.Acad.Sci.ParisSer.I, 341(2005)787–792. • K. Lipnikov, M. Shashkov, D. Svyatskiy, Yu. Vassilevski, J. Comput. Phys., 227 (2007) 492–512. The location of collocation point (cell centers) is restrictive.

8. MPFA Monotone Nonlinear monotone schemes G. W. Yuan, Z. Q. Sheng, Monotone finite volume schemes for diffusion equations on polygonal meshes, J. Comput. Phys. , 227 (2008) 6288-6312. We construct a nonlinear monotone finite volume scheme for linear diffusion problems, based on an adaptive approach of choosing stencil in the construction of discrete normal flux. The scheme is monotone on any star-shaped polygonal meshes, without any special restriction for the choice of cell centers. It is not easy to be implemented to existing code directly.

9. Repair techniques for existing schemes (2) Enforcing positivity-preserving • Priori Processing Constrained optimization R. Liska, M. Shashkov, Enforcing the discrete maximum principle for linear finite element solutions of second-order elliptic problems, Commun. Comput.Phys., 3, 4 (2008) 852–877. Nonlinear constrained finite element approximations, Limiter D. Kuzmin,M.J. Shashkov, D. Svyatskiy, Journal of Computational Physics, 228 (2009) 3448–3463.

10. (3) Enforcing negative to zero: ENZ • Posterior Processing Advantage： The common approach: “setting the negative values to be zero” is convenient to implement it in existing codes directly. Disadvantage： • Accuracy decreasing • violation of local conservation R. Liska, M. Shashkov, Commun. Comput. Phys. 3 (2008) 852-877. global conservation As we know, there exists no method of enforcing both positivity and local conservation for discrete solution of nonlinear diffusion problems.

11. Another Posterior Processingmethod (4)Returning back to previous nonlinear iteration step: Enforcing positivity algorithm (EPA) • violation of local conservation • large error Our method: returns back to previous iteration step and keeps local conservation as well.

12. (5)Take smaller time-step for nonsteady problems • too expensive computational cost • too much cumulate error Other methods: • We will not discuss the method of shortening time-step for the sake of computational cost. • In our study no more restrictive time-step is needed.

13. 2. Enforcing Positivity with Conservation Nonlinear diffusion problems

14. Picard Nonlinear Iteration The nonlinear coefficient has no definition for negative . Once the calculation is stopped.

15. Logically rectangular grid

16. Finite Volume Method

17. Nine point scheme

18. Conservative enforcing positivity algorithm：CEPA • Suppose that the solution values on all grid-cells at the s-th iteration step are nonnegative • After solving the (s+1)-th Picard iteration step, we get the solution values

19. The set of all cells is a union of three disjoint subsets: negative cells (Strictly) positive cells nonnegative cells fix up： i.e., E=P

20. then negative values are set to the values of the previous iteration step, and the fluxes are also revised to be the fluxes of previous iteration step simultaneously: ＋ ＋ Enforcing positivity for the cell inN

21. ＋ Values and flux are kept for the cell inP and the values on all its neighboring cells with common edge are also nonnegative. On (strictly) positive cells, the value and flux are kept to be the current iteration step:

22. ＋ Update the cell values in P_0 and at least one of its neighboring cells with common edge is negative cell, then at least one flux on cell-edge has been updated. On edges with neighboring cells being negative, take flux from previous iteration step; On other edges, flux remain to be current iteration step; Use the revised fluxes to renew the values on the cells in P_0: (*) It’s an explicit calculation and a key to keep local conservation.

23. Use of the five-point scheme at the first iteration step Our algorithm requires: . It implies that all the cell values must be nonnegative at least at the first iteration step. We use the classic five-point scheme (FPS) at the first iteration step to assure the numerical solutions to be nonnegative for all cells. It is easy to reduce the NPS to the FPS. We need only add a parameter in flux formula: For the first iteration step, i.e., , we choose ; otherwise .

24. The procedure of CEPA At (n+1)-th time level, algorithmic procedure can be outlined as follows: 1) Use the FPS at the first step of Picard iteration. 2) Formulate the fluxes on edges, and save them in the array FLUX_S. 3) Perform (s+1)-th iteration step and compute the fluxes, and then save them both in the array FLUX_SP1 and the array FLUX. 4) If all the cell values are positive, then go to step 7); otherwise go to next step.

25. The procedure of CEPA • 5) Scan the grid: decompose E into three disjoint subsets • For the cells with negative values (i.e., cells in N), enforce them to be equal to values of previous iteration step and replace the corresponding fluxes in FLUX with those in FLUX_S; • For the cells in P, no update is needed; • For the cells in P_0, revise their values in terms of formula (*), where is taken from FLUX, and • from FLUX_SP1;

26. The procedure of CEPA 6) Replace E with , and go to 4) Note at the step 5), update the decomposition . 7) Replace the fluxes in FLUX_S with those in FLUX, and go to 3) until the convergence criterion of nonlinear iteration is satisfied.

27. Remark1 • When we have completed the update of the solution values on E one time, in which all cell values and fluxes in N have been set to be corresponding cell values and fluxes of previous iteration step, in the following update procedure we need only to consider , and re-decompose it into . • When the updated value for becomes negative, the new N is non-empty; and we should continue to perform the above procedures until all the cells become positive cells, i.e., E=P.

28. Remark 2 • The above algorithm procedure must be stopped within finite time updates. • Notice that after performing the update procedure one time, it is unnecessary to consider the elements of the old N in all the following updates.

29. Suppose that . After updated, its surrounding cell values need to be revised for the sake of conservation. Remark 3 • Thus, the revised values may become negative, e.g., we suppose becomes negative. When we use the same way to deal with the cell , the flux between cells and is just the flux of previous iteration step, which is kept to be fixed through all the following updates. • So the cell value at will not be changed anymore.

30. Remark 4 • Therefore, we can consider only (in the step 6 of the algorithm procedure process) instead of in the following further updates. • If the new N is non-empty at this time, then we return to perform update procedures. It follows that, after updating finite times, N must become empty set and all cells are transformed to be positive cells so that E=P.

31. Remark 5 • From the construction of our algorithm of posteriori nonnegativity correction, it's easy to see that it is local conservation preserving at each step of nonlinear iteration. P N P P_N P_N P P N

32. CENZ: Conservative Enforcing Negative-values to Zero • The drawback of the standard set-to-zero fixup is that conservation is destroyed. • Another conservative fixup strategy to remedy negative values is designed. • The idea of conservative correction used in CEPA is further developed and combined with the standard set-to-zero fixup.

33. 3. Numerical Test Error Conservation error

34. Test 1：

35. x y 200200 rectangular grid

36. Nine point scheme NNPS Blankpart indicates negativecells 6464 Shestakov CEPA CENZ

37. Errors Comparison L2 conservation

38. Time : Time : x x y y CEPA NNPS

39. Test 2： Exact Solution:

40. 2020 Random Mesh

41. ENZ CENZ CPEA Error on Random Mesh

42. 4. Conclusion: Priori Processing

43. 4. Conclusion: Posterior Processing ENZ : Enforcing negative to zero EPA: Returning back to previous nonlinear iteration step CENZ:Conservative ENZ CEPA:Conservative EPA

44. 4. Conclusion • Our enforcing positivity algorithm with conservation： • positivity preserving • locally conservative, better than the common tip “setting negative values to zero” • no remarkable accuracy loss • no more restriction on time-step and mesh-cell • easy to be implemented • applicable to other finite volume schemes on structured and unstructured meshes • useful and practical in applications

45. References • Yanzhong Yao, Guangwei Yuan, Enforcing positivity with conservation for Nine-Point scheme of nonlinear diffusion equations, Comput. Methods Appl. Mech. Engrg., 2012, 223-224, 161-172. • Fujun, Cao, Yanzhong Yao, Guangwei Yuan, Positivity repair techniques with conservation for Kershaw scheme of diffusion equations, 2013.

46. Thank You！