Some improvements on Moment-of-Fluid method in 3D rectangular hexahedrons
SSome improvements on Moment-of-Fluid method in 3D rectangular hexahedrons
Zhouteng Ye a,b , Mark Sussman b , Xizeng Zhao a, ∗∗ a Ocean College, Zhejiang University, Zhoushan 316021, Zhejiang, Peoples Republic of China b Department of Applied & Computational Mathematics, Florida State University, United States
Abstract
The moment-of-fluid method (MOF) is an extension of the volume-of-fluid method with piecewise linear interfaceconstruction (VOF-PLIC). In MOF reconstruction, the optimized normal vector is determined from the referencecentroid and the volume fraction by iteration. The state-of-art work by Milcent and Lemoine [1] proposed an analyticgradient of the objective function, which greatly reduces the computational cost. In this study, we further acceleratethe MOF reconstruction algorithm by using Gauss-Newton iteration instead of Broyden-Fletcher-Goldfarb-Shanno(BFGS) iteration. We also propose an improved initial guess for MOF reconstruction, which improves the e ffi ciencyand the robustness of the MOF reconstruction algorithm. Our implementation of the code and test cases are availableon our Github repository. Keywords:
Moment of fluid, Interface reconstruction, Material tracking, Nonlinear optimization
1. Introduction
A lot of scientific and engineering problems involve with tracking the interface between di ff erent materials. Multi-ple volume tracking / capturing methods, such as volume-of-fluid method (VOF) [2], level set method [3], front trackingmethod [4] have been introduced to describe the motion of interface explicitly or implicitly. Among those methods,the volume-of-fluid method with piecewise line interface construction (VOF-PLIC) is one of the most widely usedmethods in tracking the interface within the Eulerian framework.Conventional VOF-PLIC reconstructs the normal vector of the reconstructed interface by using the a stencil thatcontains the information of the neighbouring grids, for example, Parker and Youngs’ algorithm [5], mixed Youngs-centered algorithm (MYC) [6], and the e ffi cient least squares volume-of-fluid interface reconstruction algorithm(ELVIRA) [7]. Although some of the VOF-PLIC reconstruction algorithms are second-order accuracy, when there isnot enough information from the neighbouring grids, for example, very small scale droplets, VOF-PLIC algorithmmay not reconstruct the interface accurately. ∗ Corresponding author: email: [email protected]; ∗∗ Corresponding author: email: [email protected];
Email addresses: [email protected] (Zhouteng Ye), [email protected] (Mark Sussman)
Preprint submitted to Journal of Computational Physics October 1, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p oment of Fluid method (MOF) [8, 9] introduces the centroid as the additional constraint to determine the normalvector of the reconstruction plane. When there is no data from adjacent cells used in reconstruction, MOF reconstruc-tion resolves the the interface with a smaller minimum scale than the VOF-PLIC algorithm. In MOF reconstruction,evaluating the objective function and its partial derivative is the most expensive part during iteration. The originalMOF algorithm by Dyadechko and Shashkov [8] has to call a very complex polyhedra intersection algorithm for 5times in every iteration. Although later study by Chen and Zhang [10] reduces the calling of the geometric algorithmto one time each iteration, the computational cost is still heavy. Lemoine et al. [11] made their first attempt to derive ananalytic form of that describes the objective function as the minimum distance from the reference centroid to a closed,continuous curve in 2D Cartesian grid. This is a fully 2D MOF algorithm as solution to the objective function canbe obtained by computing the cubic or quartic roots of polynomials instead of iteration. Unfortunately, this approachcannot be extended to 3D [1]. Milcent and Lemoine [1] derived an analytic form of the partial derivative of objectivefunction in 3D rectangular hexahedron, by using the analytic form, the computational cost is significantly reduced.The algorithm could be more than 200 times faster than the conventional MOF reconstruction [8].In this study, we further accelerates the MOF algorithm based on the analytic gradient by Milcent and Lemoine [1].We use the Gauss-Newton algorithm instead of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm in Milcentand Lemoine [1]. Although Gauss-Newton algorithm has been used in other MOF studies [12, 13], no detailedcomparison between the two algorithms has been carried out. We show that the Gauss-Newton algorithm significantlyreduces the number of gradient calculation. We also propose an improved form of initial guess, which provides acloser initial guess to the global minimum. The improved initial guess helps to reduce the iteration step and improvesthe robustness of the Gauss-Newton algorithm.The run-time ratio and robustness of the method could be implementation-dependent. Out implementation of thecode and test cases are available on our Github repository (https: // github.com / zhoutengye / NNMOF). All the numericaltests are done on a workstation with Intel(R) Xeon(R) Platinum 8270 processors with Intel compiler 2020.
2. Moment of Fluid reconstruction
As an extented the VOF-PLIC method, the MOF method reconstructs the interface in a 3D rectangular hextahedroncell with a plane (cid:110) x ∈ R | n · ( x − x ) + α = (cid:111) , (1)where n is the normal vector, x is the reference point of the cell, x is the origin of the Cartesian coordinate, either thecenter of the cell or the lower corner of the cell, depending on the computational algorithm. α is the parameter thatrepresents the distances from the reference point x . The volume of the cutting polyhedron by the reconstruction planesatisfies | F ref ( n , α ) − F A ( n , α ) | = . (2)2nd E MOF = (cid:107) x ref − x A ( n , b ) (cid:107) (3)In addition to the constraint on volume fraction, the MOF reconstruction also minimizes error of the centroid with E MOF ( Φ ∗ , Θ ∗ ) = (cid:107) f ( Φ ∗ , Θ ∗ ) (cid:107) = min ( Φ , Θ ):Eq . (2) holds (cid:107) f ( Φ , Θ ) (cid:107) (4)Eq. (4) minimizes the distance in 3D with two parameters by converting the normal vector in Eq. (1) to the polarangle and azimuthal angle in a spherical coordinate system. Eq. (4) is a non-linear least square problem for Φ and Θ ,it is solved with optimization algorithm via iteration. We use the Gauss-Newton algorithm to minimize Eq. (4). Starting with an initial guess of Φ , Θ , the solutionprocedure are:1. Find the centroid x k corresponds with the angle ( Φ k , Θ k ).2. Determine the Jacobian matrix J k using the analytic solution Milcent and Lemoine [1].3. Determine the sift angle ( ∆Φ k , ∆Θ k ) = − H k J Tk f k (5)where H k = J Tk J k is the Hessian matrix. In this problem, the dimension of the Hessian matrix is 2 × Φ k + , Θ k + ) = ( Φ k + ∆Φ k , Θ k + ∆Θ k ).The iteration stops while convergence conditions are full-filled. Multiple convergence criteria can be adopted,for example, centroid error, error of the gradient of the objective function, minimum advance angle, and maximumiteration step.In this problem, even though the gradient of the objective function has been significantly boosted by the analyticgradient of Milcent and Lemoine [1] compared with the conventional numerical gradient approach, calculating theobjective function f and the gradient objective function ( ∂ f ∂ Φ , ∂ f ∂ Θ ) still takes most of the computational time duringiteration. The number of calling the gradient algorithm determines the total computational cost of the iteration. Inthe original MOF method [8], the non-linear optimization Eq. (4) is solved with Broyden-Fletcher-Goldfarb-Shanno(BFGS) algorithm, which is also used in Chen and Zhang [10], Milcent and Lemoine [1]. In BFGS algorithm, theadvance angle ( ∆Φ k , ∆Θ k ) needs to be determined by a line search algorithm. In the numerical tests of Milcent andLemoine [1], every iteration needs 8-10 steps of line search, which means the total number of calling the gradientalgorithm is much bigger than the number of the iteration step.The main advantage of the BFGS algorithm over Gauss-Newton algorithm is that the BFGS algorithm approxi-mates the inverse of the Hessian matrix, which avoids the calculation of the inverse of the Hessian matrix directly.However, in this problem, the shape of the Hessian matrix is only 2 ×
2, making the cost of the inverse matrix algorithmnegligible. While in Gauss-Newton algorithm, the shift angle is directly determined by Eq. (5), so that the number of3alling the gradient algorithm equals to the iteration step. Compared with BFGS algorithm, the number of calling thegradient algorithm is much smaller than the BFGS algorithm if both algorithms converge with same iteration steps.Other non-linear optimization could potentially be used in minimizing Eq. (4). For example, the Levenberg-Manquardt algorithm [14] is known as an inproved Gauss-Newton algorithm using a trust region approach. Althoughthe Levenberg Manquardt algorithm is more robust than the Gauss-Newton algorithm, finding the value of the trustregion involves trial computation of the objective function and its gradient, which could significantly increase thecomputational cost. For e ffi ciency, we use the Gauss-Newton algorithm other than Levenberg-Manquardt algorithm.To ensure the robustness of the Gauss-Newton algorithm, we provide an improved initial guess in next subsection. The choice of initial guess is important because there may exists multiple local minima in the objective function.Dyadechko and Shashkov [8] suggested the form n = x c ( Ω ) − x ref (6)as safe initial guess. Dyadechko and Shashkov [8] also claimed that the line-search algorithm guarantees the initialguess finally reaching the global minima. In Gauss-Newton iteration, the step is automatically determined, there is notrial-step selection. The Gauss-Newton algorithm is more likely to be sensitive to the initial guess than BFGS withline search algorithm used in the study of Dyadechko and Shashkov [8], Milcent and Lemoine [1]. (a) Evolution with fixed angle . . . . . . . . . . . . (b) Locus of centroids (c) Evolution with fixed volume fractionFigure 1: The locus of the centroids for reference with fixed volume fraction V (Dashed line)and the locus of cenroids with fixed azimuthal angle Φ (Solid lines) in a unit 2D cell.The black curve corresponds with the locus of centroids for reference volume V = . We propose a new form of the initial guess in this section. To better demonstrate the philosophy of our proposedinitial guess, we simplify the 3D problem to 2D by setting the polar angles
Φ = π/
2, which simplifies the 3D problemto a 2D problem. Fig. 1 shows the locus of the centroids of the cutting polygon by a line interface in a unit cell. Thesolid lines in Fig. 1 (b) corresponds with the evolution of the centroid with a fixed azimuthal angle Φ (See Fig. 1(a)),4nd the dashed lines in Fig. 1(b) corresponds with the evolution of the centroid with a fixed volume fraction V (SeeFig. 1(c)). When the volume fraction V > / V > / x of an initial guess Φ could bevery close to the reference centroid x c , but has a big error with the exact azimuthal angle Φ .We show the error of initial guess Eq. (6) in Fig. 2 (a). Eq. (6) gives a good initial guess in most areas exceptfor the region that is near the face of the cell. Those regions correspond with very small volume fraction. We proposeanother candidate initial guess by assuming the reference centroid as the centroid of a right triangle reconstruction (ora trirectangular tetrahedron in 3D). n = v ( Ω ) − x ref , (7)where ˜x v ( Ω ) is the vertex of the quadrant (or octant in 3D) that x ref is located. The error map of the azimuthal angleis plotted in Fig. 2 (b). The right triangle approximates the small volume correctly especially when the centroid isnear one of the vortex of the grid cell. We evaluates the value of the objective function with the two candidate initialguesses and pick the one with smaller error of the centroted n = arg min { n , n } (cid:110) E MOF (cid:16) n (cid:17) , E MOF (cid:16) n (cid:17)(cid:111) (8)The error of the azimuthal angle error ∆Φ e of Eq. (8) are plotted in Fig. 2 (c). In 2D case, the maximum error of thepolar angle by Eq. (8) is approximately π/
20, while the maximum error of the polar angle from Eq. (6) is about π . Wealso teste the initial guess in 3D. The error of the initial guess ∆Θ + ∆Φ by Eq. (6) is about π . By using our improvedinitial guess, the error reduces to about π . . . . . . . . . . . . . . . . . . (a) Initial guess Error of Eq. (6) . . . . . . . . . . . . . . . . . (b) Initial guess Error of Eq. (7) . . . . . . . . . . . . . . . . . (c) Initial guess Error of Eq. (8)Figure 2: The azimuthal angle error ∆Θ between the initial guess and exact angle. . . . . . . C C o un t o f d a t a s e t (a) Exponential case . . . . . . C C o un t o f d a t a s e t (b) Uniform case . . . . . . C C o un t o f d a t a s e t (c) Extreme caseFigure 3: Distributions of volume fraction for three data sets.Table 1: tab: L error of the reconstruction test Gauss-Newton BFGSexponential uniform extreme exponential uniform extremeOriginal initial guess 1.98e-09 2.92e-06 1.19e-05 2.08e-09 2.27e-09 2.74e-09Improved initial guess 1.97e-09 1.93e-09 1.85e-09 2.08e-09 2.01e-09 1.93e-09
3. Numerical tests
In this section, we test the accuracy and robustness of our MOF reconstruction with Gauss-Newton algorithm withimproved initial guess. Two criteria are evaluated: the CPU time and the robustness. Three data sets are generated byfinding the exact centroid of a cutting polyhedron of a unit cube by a plane. We use data-sets with di ff erent distributionto show the performance of our algorithm, especially the robustness for extreme cases (See Fig. 3): (1) Exponentialcase with a normal distribution; (2) Uniform case corresponds with uniform distribution; (3) Extreme case with ashifted normal distribution which contains more values near 0 or 1. In this test, the tolerance for iteration is 10 − , themaximum iteration step is 100.In Table 1, the error of the Gauss-Newton algorithm with original initial guess increases when more extremedata appears in the test data set, while the BFGS shows a better robustness than the Gauss-Newton algorithm. Withthe improved initial guess, both algorithms show a very good robustness in all test cases. In BFGS algorithm, each Table 2: Run-time ratio and iteration steps in reconstruction test
Gauss-Newton BFGSAveragediteration Averagedline search run-timeratio Averagediteration Averagedline search run-timeratioOriginal initial guess 3.34 0 1.38 6.26 12.56 5.29Improved initial guess 2.48 0 1 5.4 10.05 4.256 able 3: Run time ration and iteration in advection test (Zalesak’s problem)
Gauss-Newton BFGSAveragediteration Averagedline search run-timeratio Averagediteration Averagedline search run-timeratioGrid: 50 × × × × In the previous test, the optimized centroid is consistent with the reference centroid. However, in practical, theoptimized centroid may not be consistent with the reference centroid. In order to test the accuracy and robustness ofthe proposed algorithm with non-linear reconstruction, we test our algorithm with a 3D Zalesak’s problem Enrightet al. [15]. We use a directional splitting Lagrangian Explicit scheme for the advection of volume fraction [6] andupdates the centroid by calculating the evolution of corresponding Lagrangian centroid. For advection of the volumefraction and centroid, We use a directional splitting Lagrangian Explicit scheme similar with the VOF-PLIC advectionin Aulisa et al. [6].The di ff erence between di ff erent method are very small (With an L error of O − ), which shows the robustness ofour algorithm. The averaged time of iteration and computational cost are listed in Table 3. With our improved initialguess, the averaged number of iteration decreases in Gauss-Newton algorithm and BFGS algorithm. The Gauss-Newton algorithm with the improved initial guess gets an acceleration of 3 to the algorithm of the [1].7 . Conclusions In this study, we show that using Gauss-Newton algorithm instead of BFGS algorithm significantly helps to accel-erate the iteration in MOF reconstruction. We also proposed an improved initial guess which makes the Gauss-Newtoniteration more robust. Our improved initial guess along with the Gauss-Newton algorithm is about 4 times faster thanthe BFGS algorithm by Milcent and Lemoine [1] in reconstruction and about 2 times faster in advection test.
5. Acknowledgments
The support provided by National Science Foundation of China (Grant Nos. 51979245, 51679212) and ChinaScholarship Council (CSC) and the during a visit of Zhouteng Ye to Florida State University is acknowledged.
References [1] T. Milcent, A. Lemoine, Moment-of-fluid analytic reconstruction on 3D rectangular hexahedrons, Journal of Computational Physics 409(2020) 109346.[2] C. Hirt, B. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1981)201–225. Number: 1.[3] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, Journal ofcomputational physics 79 (1988) 12–49. Number: 1.[4] S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of Computational Physics 100(1992) 25–37.[5] B. Parker, D. Youngs, Two and three dimensional Eulerian simulation of fluid flow with material interfaces, Atomic Weapons Establishment,1992.[6] E. Aulisa, S. Manservisi, R. Scardovelli, S. Zaleski, Interface reconstruction with least-squares fit and split advection in three-dimensionalCartesian geometry, Journal of Computational Physics 225 (2007) 2301–2319. Number: 2.[7] J. E. Pilliod, E. G. Puckett, Second-order accurate volume-of-fluid algorithms for tracking material interfaces, Journal of ComputationalPhysics 199 (2004) 465–502. Number: 2.[8] V. Dyadechko, M. Shashkov, Moment-of-fluid interface reconstruction, Los Alamos Report LA-UR-05-7571 (2005).[9] V. Dyadechko, M. Shashkov, Reconstruction of multi-material interfaces from moment data, Journal of Computational Physics 227 (2008)5361–5384.[10] X. Chen, X. Zhang, An improved 3D MoF method based on analytical partial derivatives, Journal of Computational Physics 326 (2016)156–170.[11] A. Lemoine, S. Glockner, J. Breil, Moment-of-fluid analytic reconstruction on 2D Cartesian grids, Journal of Computational Physics 328(2017) 131–139.[12] M. Jemison, E. Loch, M. Sussman, M. Shashkov, M. Arienti, M. Ohta, Y. Wang, A Coupled Level Set-Moment of Fluid Method forIncompressible Two-Phase Flows, Journal of Scientific Computing 54 (2013) 454–491. Number: 2-3.[13] A. Asuri Mukundan, T. Mnard, J. C. Brndle de Motta, A. Berlemont, A 3D Moment of Fluid method for simulating complex turbulentmultiphase flows, Computers & Fluids 198 (2020) 104364.[14] K. Madsen, H. B. Nielsen, O. Tingle ff , Methods for non-linear least squares problems (2004).[15] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A Hybrid Particle Level Set Method for Improved Interface Capturing, Journal of Computa-tional Physics 183 (2002) 83–116. Number: 1., Methods for non-linear least squares problems (2004).[15] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A Hybrid Particle Level Set Method for Improved Interface Capturing, Journal of Computa-tional Physics 183 (2002) 83–116. Number: 1.