High-Order Kinetic Relaxation Schemes as High-Accuracy Poisson Solvers
aa r X i v : . [ phy s i c s . c o m p - ph ] A p r November 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC
High-Order Kinetic Relaxation Schemesas High-Accuracy Poisson Solvers
M. Mendoza
ETH Z¨urich, Computational Physics for Engineering Materials, Institute for BuildingMaterials, Wolfgang-Pauli-Str. 27, HIT, CH-8093 Z¨urich, [email protected]
S. Succi
Istituto per le Applicazioni del Calcolo C.N.R., Via dei Taurini, 19 00185, Rome, [email protected]
H. J. Herrmann
ETH Z¨urich, Computational Physics for Engineering Materials, Institute for BuildingMaterials, Schafmattstrasse 6, HIF, CH-8093 Z¨urich, [email protected]
We present a new approach to find accurate solutions to the Poisson equation, as ob-tained from the steady-state limit of a diffusion equation with strong source terms. Forthis purpose, we start from Boltzmann’s kinetic theory and investigate the influence ofhigher order terms on the resulting macroscopic equations. By performing an appropriateexpansion of the equilibrium distribution, we provide a method to remove the unneces-sary terms up to a desired order and show that it is possible to find, with high level ofaccuracy, the steady-state solution of the diffusion equation for sizeable Knudsen num-bers. In order to test our kinetic approach, we discretise the Boltzmann equation andsolve the Poisson equation, spending up to six order of magnitude less computationaltime for a given precision than standard lattice Boltzmann methods.
Keywords : High Knudsen number, higher-order moments, diffusion equation, Poissonequation, lattice Boltzmann
1. Introduction
The diffusion equation is widely used to describe transport phenomena in whichheat, mass, and other physical quantities are transferred in space and time due tounderlying molecular collision processes 1. At steady-state and when the sourceterm does not depend explicitly on the concentration of the associated field, suchdiffusion processes settle down in the form of a non-local relation between the spatialdistribution of the source and the resulting concentration of the associated field, arelation which is governed by the Poisson equation 2. The notion of the Poissonequation as the steady-state solution of the diffusion equation extends to many otherphenomena, e.g. electrostatics systems 3, plasma physics 4, chemistry and biology5 ,
6, molecular biophysics 7, density functional theory and solid-state physics 8 , , ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC to name but a few. In fact, solving the Poisson equation accurately and efficientlycontinues to be a subject of intense research to this day.From a microscopic point of view, the diffusion equation emerges from the un-derlying Boltzmann kinetic equation in the limit of very small Knudsen numbers,i.e. molecular mean free path much shorter than the typical macroscopic scale. Solv-ing the Boltzmann kinetic equation to study diffusion processes makes apparentlylittle sense, since the Boltzmann distribution function lives in a double-dimensionalphase-space, i.e. position and velocity. However, in the last two decades, minimal(lattice) versions of the Boltzmann equation have been developed, in which veloc-ity space is reduced to a small set of discrete velocities, so that the computationalcost is cut down dramatically, making the solution of the lattice kinetic equationoften competitive towards standard grid-discretisation of the corresponding partialdifferential equation.The Lattice Boltzmann Method (LBM) has attracted much interest for solvingthe Navier-Stokes equations 11 ,
12, and has been extended to solve other type ofphenomena, e.g. diffusion equation 13 ,
14, Maxwell equations15, quantum systems16, relativistic fluid dynamics 17 ,
18, and the Poisson equation 2 , , , ,
22, to namebut a few. To the best of our knowledge, to date, the solution of the Poisson equationbased on kinetic theory has been confined to smooth (non-stiff) source terms andsmall Knudsen numbers 23. If the target is the steady-state solution of the diffusionequation, as it is the case for the Poisson equation, the kinetic approach loosesmuch of its appeal, since its inherently real-time dynamics must proceed throughsmaller time-steps than those affordable by fictitious-time iteration techniques. Inthis paper we shall show that such gap can be closed by resorting to higher-orderkinetic formulations, at least for the cases where the source term does not dependon the concentration field avoiding artificial solutions as pointed out in Ref. 2.Higher-order formulations of the LB equation have been developed before,mostly in the context of thermal 24 and relativistic fluid dynamics 25, and recently,density functional theory 26. There, the main idea is to design the equilibrium dis-tribution in such a way to include not only the conserved moments, such as densityand current, but also the non-conserved ones, namely the momentum flux tensorand the heat flux, and eventually higher order kinetic moments with no direct hy-drodynamic significance (sometimes called “ghosts”). Clearly, in order to matchthis longer ladder of kinetic moments, a correspondingly larger set of discrete ve-locities is required. It should be emphasized that the discrete kinetic equation is asuperset of the diffusion equation, and consequently, to the purpose of solving thediffusion equation alone, higher-order contributions must be regarded as discretisa-tion artefacts which need to be minimised, and possibly canceled altogether. This isprecisely the target of this paper. Here, we will use higher-order lattices and associ-ated equilibria, to remove higher-order contributions to the macroscopic equationsreproduced by the Boltzmann equation. For this purpose, we develop the respectiveLBM and show that, by expanding the equilibrium distribution up to seventh orderin Hermite polynomials, one can solve the Poisson equation six orders of magnitudeovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC faster than with standard LB models, at a given level accuracy.The paper is divided as follows: In Sec. 2, we introduce the diffusion equationwith stiff source term and the associated Boltzmann kinetic equation. In Sec. 3, wedescribe the influence of higher order moments on the steady-state solution of thediffusion equation. In Sec. 4, we develop a LBM for the Poisson equation, validateand compare the model with standard cases. Finally, we discuss the results andoutlook of our work.
2. From Boltzmann kinetic theory to diffusion equation
Let us begin by writing the diffusion equation for a scalar Φ in the standard form: ∂ Φ ∂t = D ∇ Φ +
S , (1)where Φ is the concentration of particles (temperature in the case of the heat trans-fer equation), D is the diffusivity, and S is the source term, which, in our case,depending neither on Φ nor on time. The steady-state version of Eq. (1) deliversthe Poisson equation with source − S/D . It is well-known that this equation can beobtained from an underlying Boltzmann kinetic equation, via a Chapman-Enskogexpansion, in the limit of small Knudsen numbers, Kn = λ/L , being λ the mean-freepath and L the characteristic size of the system 27.The Boltzmann equation in the single relaxation time approximation (BGK) 28reads as follows: ∂f∂t + ~v · ∇ f = − τ ( f − f eq ) + S , (2)where f ≡ f ( ~x, ~v, t ) is the probability distribution function, ~v are the microscopicvelocity vectors, τ is the relaxation time, and f eq is the equilibrium distribution.In the non-relativistic context, the local equilibrium is given by a Maxwell-Boltzmann distribution: f eq = e − ( v − u ) / c s (2 πc s ) / , (3)where u is the local flow speed. Such local equilibrium encodes the basic symmetriesof Newtonian mechanics, particularly Galilean invariance, which is built-in via thedependence on the molecular speed relative to the fluid v − u , rather than theabsolute one, v . In order to secure Galilean invariance for any fluid speed, kineticmoments at all orders should be matched on the lattice. This would require aninfinite connectivity, which is clearly unviable on any realistic lattice. As a result,LB local equilibria are typically based on finite-order Hermite expansions, of theform 29: f eq ( ~x, ~v, t ) = w ( ~v ) ∞ X n =0 a ( n ) n ( ~x, t ) H ( n ) n ( ~v ) , (4)ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC where H ( n ) n are the tensorial Hermite polynomials of order n , and w ( v ) is the weightfunction defined as, w ( ~v ) = 1(2 πc s ) / e − v / c s , (5) c s being the sound speed. Note that w ( v ) is the uniform global equilibrium, corre-sponding the no-flow limit of the local equilibrium f eq . The coefficients a ( n ) n (kineticmoments) are calculated by projecting the equilibrium distribution onto the respec-tive Hermite polynomial, a ( n ) n = Z f eq H ( n ) n d v . (6)The lowest-order kinetic moments have a direct macroscopic interpretation. Forinstance, by integrating on the velocity space ( H = 1),Φ = Z f d v = Z f eq d v , (7)where the last equality comes from the requirement of mass conservation on thecollsion operator, namely: 1 τ Z ( f − f eq ) d v = 0 . (8)The diffusivity of the system is dictated by the relaxation time as D = c s τ. (9)It can be shown that the Boltzmann equation recovers the diffusion equationin the limit of small Knudsen numbers, which is related with the mean-free pathand therefore with the relaxation time. This is the strongly-interacting fluid regime,whereby the particle dynamics is collision-dominated, so that collective behavioursets in on a short timescale τ . The opposite limit of large Knudsen numbers, hencelarge values of τ , corresponds to the free-particle motion (ballistic regime), wherebymemory of the initial conditions is kept for a very long time. It is therefore clear thatin the ballistic regime higher order moments of the equilibrium distribution do notrelax and consequently the Boltzmann equation does not converge to any standarddiffusion equation. On the other hand, by increasing the relaxation time, hence theeffective diffusivity, one would intuitively expect a quickest path to steady-state. Inthis respect, it would be highly desirable to derive a diffusion equation in kineticform, which can attain steady-state solution at relatively large Knudsen numbersby exactly cancelling as many high-order terms as possible.As discussed above, this can be achieved by zeroing as many higher order con-tributions as possible, while retaining the largest possible diffusion coefficient.ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC
3. Influence of Higher-Order Terms
Let us consider the case when all time derivatives are zero, ∂/∂t = 0. The Boltzmannequation, Eq. (2), can be written as, f = f eq + τ S − τ~v · ∇ f . (10)By integrating this equation in velocity space, we obtain a steady-state continuityequation for the current density: ∇ · ~J = S , (11)where ~J = R f~v d v and S = R S d v .Multiplying Eq. (10) by ~v and integrating again, we obtain ~J = ~J eq + τ ~S − τ ∇ · Π , (12)with Π ≡ Π αβ = R f v α v β d v being the momentum-flux tensor and ~S = R S ~v d v the source current. The first term at the right-hand-side takes the form ~J eq = Φ ~u ,where ~u is the flow speed in the local equilibrium. For the case of a purely diffusivedynamics this term is zero. The second term contributes a shift τ ~S to the flow speed,and it must also be set to zero for the case of pure diffusion. Finally, the pressuretensor shoud reduce to c s Φ, so that inserting (12) into (11), the diffusion equationis obtained. The first two conditions are automatically ensured by setting ~u = 0 inthe local equilibrium. The third one, however, cannot be enforced exactly becauseof higher order contributions to the momemtum flux tensors.Indeed, by iterating the procedure, one obtains the general expression: ∞ X n =0 ( − n τ n ∇ n +1 h Π eq( n +1) + τ S ( n +1) i = S , (13)where Π eq( n ) ≡ R f eq v α ...v α n d v and S ( n ) ≡ R S v α ...v α n d v , are both tensorsof order n . Note that for small Knudsen numbers ( τ ∇ ≪
1) higher-order termsbecome negligible and the series is convergent. We also observe that the solution isdictated by the equilibrium distribution and the source term, so that, by properlychoosing both expressions, higher order terms can be canceled thus opening the wayto high-accuracy solutions of the Poisson equation.In numerical practice, accuracy is typically improved by increasing the resolu-tion of the grid. This imposes correspondingly smaller time steps, which is clearlyexpensive, especially in three dimensions. High-order methods are meant to mit-igate the problem by achieving the same level of accuracy with many less gridpoints. They are typically based on higher-order stencils for the discretisation ofthe corresponding differential operators.In this work, we present a procedure whereby the same goal is achieved bydrawing directly from an underlying lattice kinetic theory, i.e. by a suitable (non-Gaussian) extension of the local kinetic equilibria as combined with the use of largersets of discrete speeds than those typically employed in standard Lattice Boltzmanntheory.ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC
4. Kinetic approach to the Poisson equation
In order for the Boltzmann kinetics to converge to the Poisson equation, the follow-ing conditions need to be met:Π eq(1) = Π eq( n> = 0 , (14) S ( n> = 0 . The equilibrium distribution can be expressed as a separable product of a functionof the microscopic velocity (global uniform equilibrium) and a function of the time-spatial coordinates, namely: f eq = φ ( ~v )Φ( ~x, t ) , (15)and similarly for the source term: S = χ ( ~v ) S ( ~x, t ) . (16)In the above, χ ( ~v ) and φ ( ~v ) are velocity dependent functions that need to be de-termined based on the conditions (14) above. To this purpose, we write: φ ( ~v ) = w ( ~v ) N X n =0 b ( n ) n H ( n ) n ( ~v ) , (17)and, χ ( ~v ) = w ( ~v ) N X n =0 c ( n ) n H ( n ) n ( ~v ) , (18)where b ( n ) n and c ( n ) n are constant coefficients and N is a cut-off truncating the Her-mite expansion. With reference to the one-dimensional case and an expansion upto seventh order ( N = 7), we obtain: φ ( v ) ≃ w ( v ) (cid:18) − c s − c s v + v c s − c s − c s v + 15 c s v − v c s (cid:19) , (19)and χ ( v ) ≃ w ( v ) (cid:18) (cid:20) − v c s (cid:21) + 3 c s − c s v + v c s + 15 c s − c s v + 15 c s v − v c s (cid:19) , (20)where the weight w ( v ) for the one-dimensional case is defined as w ( v ) = 1(2 πc s ) / e − v / c s . (21)With Eqs. (19) and (20), one can prove that the the following expressions for themoments are fulfilled: Π eq(0) = Φ, Π eq(1) = 0, Π eq(2) = Φ c s , Π eq( n> = 0, S (0) = S , S ( n> = 0. These expressions hold only up to N = 7.To develop a lattice Boltzmann model and solve the Poisson equation numeri-cally, we need to impose a quadrature and find a corresponding set of finite velocityovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC vectors, v i , such that the orthogonality conditions are fulfilled to the desired order24 ,
25, namely: Z w ( v ) H n ( v ) H m ( v ) dv = N v X i =0 w i H n ( v i ) H m ( v i ) . (22)Here w i are discrete weights and H n are the one-dimensional Hermite polynomials.The number of velocity vectors N v depends on the order of the expansion that weintend to reproduce with the quadrature. For expansions up to seventh order, weneed at least 13 velocity vectors and weights. Each set of numbers also provides itsown lattice sound speed c s , which goes from about 0 . N = 3 to nearly 1 . N = 7. The details of these values are given in Appendix A.Higher-order schemes are usually exposed to numerical instabilities, due to theappearance of additional modes in the dispersion relation. However, since our pro-cedure is designed to annihilate precisely these modes, we expect our approach tobe stable. Furthermore, the method proposed here can also be used as a systematictechnique of calculating weights and velocity vectors, such that one can achieve thedesired level of accuracy for the computation of Laplacian operators, similar to thework proposed in Refs. 30 , φ ( ~v ) ≃ w ( ~v ) (cid:18) − c s − c s ~v + ~v c s (cid:19) , (23)and χ ( ~v ) ≃ w ( ~v ) (cid:18) (cid:20) − ~v c s (cid:21) + 15 c s − c s ~v + ~v c s (cid:19) , (24)where the weight w ( ~v ) is defined as w ( ~v ) = 1(2 πc s ) / e − ~v / c s . (25)For the three-dimensional case, we should calculate the corresponding set of discretevelocity vectors, ~v i , such that the orthogonality condition is fulfilled up to fifth order,namely: Z w ( ~v ) H ( m ) m ( ~v ) H ( n ) n ( ~v ) d v = N v X i =0 w i H ( m ) m ( ~v i ) H ( n ) n ( ~v i ) , (26)which is the three-dimensional version of Eq. (22). By solving these algebraic equa-tions, we find that we need at least 111 velocity vectors and 10 weights. Details aregiven in Appendix A.ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC -500-400-300-200-100 0 100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 φ x/L N = 1N = 3N = 5N = 7Theory -0.1-0.08-0.06-0.04-0.02 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 ∆ φ x/L Fig. 1. Solution of the Poisson equation for different cut-off N and l = 1. The absolute error iscomputed ∆ φ ≡ φ n − φ t , where φ n is the potential calculated with the lattice Boltzmann model,and φ t is the analytical solution of Eq. (27). The system size is L = 128 lattice cells.
5. Numerical Results
As an application, we solve the Poisson equation 3, ∇ Φ = − ρǫ , (27)where ρ is the charge density and ǫ is the electric permittivity.Replacing the equilibrium distribution f eq from Eq. (15) into Eq. (13), we obtain: S + c s τ ∇ Φ ≃ , (28)implying that S = ρc s τ /ǫ . In the lattice Boltzmann model, one must take intoaccount second order corrections to the relaxation time, τ = τ lb − /
2, where τ lb thelattice Boltzmann relaxation time.From Eq. (13), it is appreciated that the spatial derivatives of the source alsointroduce errors in the solution. Thus, by annihilating higher order contributions,we are also eliminating spurious effects due to the derivatives of the source term. Ifthe source term is smooth, this derivatives are negligible, however, for stiffer sourceterms, one can use the present approach to add accuracy to the solution. In orderto investigate this issue, we solve the potential for the following charge density: ρ ( x ) = ǫ sin (cid:18) πxL l (cid:19) , (29)where the integer l controls the smoothness of the derivatives of the charge density,and L is the simulated length, with x ∈ [0 , L ). We use periodic boundary conditionsfor simplicity, and τ lb = 1 (numerical units).We express the relative error as ǫ = aL − α (30)where the amplitude a and scaling exponent α (order of accuracy) both depend on N . From Fig. 1, we see that for a relatively large system size, L = 128, the first orderexpansion is sufficient to produce satisfactory results, with errors around 0 .
02% forovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC -30-20-10 0 10 20 30 40 0 0.2 0.4 0.6 0.8 1 φ x/LN = 1N = 3N = 5 N = 7Theory -0.1-0.08-0.06-0.04-0.02 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 ∆ φ x/L Fig. 2. Solution of the Poisson equation for different cut-off N and l = 4. The absolute error iscomputed ∆ φ ≡ φ n − φ t , where φ n is the potential calculated with the lattice Boltzmann model,and φ t is the analytical solution of Eq. (27). The system size is L = 128 lattice cells. -2-1 0 1 2 3 0 0.2 0.4 0.6 0.8 1 φ x/LN = 1N = 3N = 5 N = 7Theory -0.1-0.08-0.06-0.04-0.02 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 ∆ φ x/L Fig. 3. Solution of the Poisson equation for different cut-off N and l = 4. The absolute error iscomputed ∆ φ ≡ φ n − φ t , where φ n is the potential calculated with the lattice Boltzmann model,and φ t is the analytical solution of Eq. (27). The system size is L = 32 lattice cells. N = 1 and N = 3, and below 10 − % for N = 5 and N = 7. In Fig. 2, we see thatby increasing the derivatives of the charge density (increasing l ), the error increasesfor all values of the cut-off N . However, as we have mentioned before, the effect ofremoving higher-order errors becomes crucial for the case of small system sizes. Bydecreasing the number of lattice cells (see Fig. 3), the discrepancies become visualand the errors are around 5% for N = 1 ,
3, 1% for N = 5, and less than 0 .
4% for N = 7.Finally, by further reducing the lattice resolution, L = 8 cells, and consideringonly one oscillation, l = 1, we see that the errors for the first order expansion, N = 1, are around 7%, while for the highest values of N = 7 they remain below0 . -2-1.5-1-0.5 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 φ x/L N = 1N = 3N = 5N = 7Theory -0.1-0.08-0.06-0.04-0.02 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 ∆ φ x/L Fig. 4. Solution of the Poisson equation for different cut-off N and l = 1. The absolute error iscomputed ∆ φ ≡ φ n − φ t , where φ n is the potential calculated with the lattice Boltzmann model,and φ t is the analytical solution of Eq. (27). The system size is L = 8 lattice cells. −9 −8 −7 −6 −5 −4 −3 −2 −1 L (cells) ε N = 1N = 3N = 5N = 7 α = 1 α = 3 α = 5 L (cells) M N = 1N = 3N = 5N = 7
Fig. 5. (left) Relative error as a function of the system size L for each cut-off value N of theHermite polynomials expansion. The lines denote the fitting curves using the general expression, ǫ = a N L − α (right). Number of iterations as a function of the system size L for each cut-off valueof N . The lines denote the fitting curve using the general expression, M = b N L , M being thenumber of time-steps. higher order lattices also bring aditional side-benefits such as a higher sound speed c s , leading to a larger diffusion coefficient and, consequently, to a faster path to thesteady-state solution. In order to highlight the concrete advantages of our approach,we next implement the same simulations for different orders N and inspect the CPUtime, number of iterations, and system size, required to achieve a specific level ofaccuracy. Computational performance
Let us next change the order of the expansion and analyse how the relative errordecreases with the system size. From Fig. 5 (left), we see that, according to ourexpectations, by increasing the order of the expansion, the order of convergence,ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC a N and b N for each value of N . Alsoshown are the accuracy exponent α N and the number oflattice sites updated per second, ξ N , in Msites/s (Millionsites per second).N a N b N ξ N (Msites/s) ξ N /b N α N .
645 10 .
98 26 . .
39 13 1 .
645 10 .
98 21 . .
95 15 33 .
11 5 .
46 5 .
44 0 .
997 37 992 . .
23 5 .
76 1 .
78 5 measured by the exponent α , also increases and in a very substantial way. Fromthe same figure (right), we can also observe that for all values of N , the numberof time-steps M scales quadratically with the lattice size, and decreases with theorder of the expansion N . The former feature is the typical signature of diffusionprocess, while the second reflects the fact that higher order lattice deliver a largersound speed, hence they take less time-steps to achieve a given level of accuracy.The computational time required to achieve a given accuracy clearly depends onthe computational speed of the model on each given machine. This is convenientlyexpressed in terms of the number of lattice sites updates per CPU second, ξ N . Thisquantity is expected to be a decreasing function of N , since higher-order quadraturesrequire more operations per lattice site.Actual measurements give the values reported in Table 1. Although no cleartrend can be extracted from these measurements, it is observed that the cases N = 5 and N = 7 are about five times more expensive than the cases N = 1 and N = 3. Also to be noted that ξ > ξ , which is counter-intuitive. However, uponinspecting the velocity vectors (see Appendix Appendix A), one can appreciate thatfor N = 7 the velocities are consecutive (no velocity gap), while for N = 5 theyare not. Since, memory access is faster for consecutive elements, we expect the case N = 7 to perform better than N = 5. More generally, the dependence on N appearsto be highly dependent on memory access patterns, which are not easily quantifiedthrough simple scaling relations.The total CPU time, t L,N , it takes to attain steady-state for a given system size L can be written as follows: t L,N = b N ξ N L . (31)In the above, the cubic dependence derives from the L diffusive scaling, as combinedwith the number of sites L (the space-time computational volume of a diffusiveprocess scales like L D in D spatial dimensions). Finally, b N is a relative scale forthe number of time-steps, which is expected to decrease at increasing N becauseof the increasing sound speed. The actual values of obtained in the simulationsare reported in Table 1. As an example, for the case L = 100, we find that thesimulations take 0 .
42 (65629 time steps), 0 .
51 (65629 time steps), 1 . .
56 (20079 time steps) seconds, for N = 1 , , , −4 −2 N t ( s ) , L ( ce ll s ) t(s)L (cells) Fig. 6. Computational time t N and system size L needed to recover the solution of the Poissonequation with a relative error of 0 . N . all values of N , while the error with respect to the analytical solution decreasesdramatically at increasing N (see Fig. 5). Hence, we can conclude that increasing N leads to a dramatic boost of accuracy at a very moderate extra computationalcost. Finally, we have inspected the computational time to steady-state for a givenaccuracy, at changing the size of the problem and the order of the quadrature. Bycombining the previous relations, (30) and (31), we obtain: t N ( ǫ ) = b N ξ N (cid:16) a N ǫ (cid:17) /α N . (32)The plot reports the computational time and system size needed to achieve a pre-scribed relative error, for the case in point ǫ = 10 − . From fig. 6, we see that byincreasing the order of the expansion in Hermite polynomials, it is possible to reducethe computational time by several orders of magnitude (up to six).These data refer specifically to the solution of the Poisson equation, but we havereasons to believe that similar conclusions would hold for other partial differentialequations with a kinetic-theory background, including many linear non-linear partialdifferential equations of great relevance in many branches of modern science.5.1.1. Comparison with multi-grid methods
In order to gain information about the performance of this approach, as comparedwith existing tools to solve the Poisson equation, we have implemented simulationsusing the multi-grid (MG) method 32 ,
33. The number of nodes in the MG simu-lations is denoted by L , like for the lattice Boltzmann method. Fixing l = 1 (fromEq. (29)), in Table 2 we show the relative errors per node, ǫ/L , together with thecorresponding computational time. Note that for N = 3, the MG method is fasterthan our approach, with the same second order accuracy. By increasing N , however,the MG becomes increasingly faster, but also increasingly less accurate at a givengrid size L . For instance, in order to attain the same accuracy of the case N = 7 onovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC ǫ/L , for the solution of the Pois-son equation using different orders N and the multi-grid method ( MG ).Inside the parenthesis, we report the computational time in milliseconds. L N = 3 N = 5 N = 7 MG
16 10 − (2) 9 × − (2) 10 − (2) 10 − ( < × − (4) 6 × − (4) 2 × − (6) 3 × − ( < × − (26) 4 × − (29) 3 × − (38) 8 × − ( < − − − × − (12) L = 64, the MG solver requires L = 16384 grid points, being still faster by abouta factor 3. Given that N = 7 requires 13 arrays, the LB memory saving is abouta factor 20. These figures indicate that the present high-order LB Poisson solver isconsistently slower than MG, however it also takes correspondingly less memory toachieve high levels of accuracy. Given that MG represents the state-of-the art forfast Poisson solvers, these can be regarded as satisfactory results.In summary, we have shown that the present LB approach provides a viableoption for the solution of Poisson equation, in the sense of simplicity and possiblyalso in terms of memory demand at high accuracy.
6. Conclusions
Summarizing, we have studied the effect of higher order terms on the steady-statesolution of the diffusion equation, using a lattice version of the Boltzmann equa-tion with high-order (non-gaussian) equilibria. For the numerical solution, we havedeveloped a higher order lattice Boltzmann model capable of reproducing up toseventh order moment of the equilibrium distribution and source charge. This per-mits to cancel exactly error terms up to seventh order. For the validation and studyof our approach, we have solved the one-dimensional Poisson equation, and foundthat the computational time to steady-state, at a prescribed level of accuracy, canbe reduced by up to six orders of magnitude as compared to standard LB formu-lations. The six orders gain in computational time is due to the fact that higherorder lattices permit to utilize much smaller grid sizes. Moreover, since they containa large number of discrete velocities, they allow faster propagation within the lat-tice, hence a correspondingly faster attainement of the steady-state. We have alsocompared with state-of-the art multigrid solvers, and found that high-order LB isnearly competitive in efficiency at a lower memory demand.Extensions to two and three-dimensions, as well as time-dependent diffusion andadvection-diffusion equations will make the object of future research.
Acknowledgements
We acknowledge financial support from the European Research Council (ERC) Ad-vanced Grant 319968-FlowCCS.ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC Appendix A. Quadrature of the Boltzmann Equation
Using Eq. (22) for each order of expansion N , we can calculate the velocity vectors v i with the respective weights w i . In addition, each set of w i and v i provides acharacteristic normalised speed c s , c s = P Mi =0 w i v i . Here we only show the casesfor N = 3 ,
5, and 7, since the simulations for N = 1 were performed with the samelattice than for N = 3. For N = 3, we have w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . v = 0 , v = 1 , v = − , v = 3 , v = − c s = 1 − p / , (A.1)for N = 5, w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 8 . × − w = 8 . × − v = 0 , v = 1 , v = − , v = 2 , v = − ,v = 3 , v = − , v = 5 , v = − c s = 0 . , (A.2)ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC ~v i with their corresponding weights w i for N = 5. There are 111 velocity vectors and 10 different weights. w i ~v i . , , . ± , , , (0 , ± , , (0 , , ± . ± , ± , , (0 , ± , ± , ( ± , , ± . × − ( ± , ± , ± . × − ( ± , ± , , ( ± , , ± , (0 , ± , ± . × − ( ± , ± , , ( ± , , ± , (0 , ± , ± ± , ± , , ( ± , , ± , (0 , ± , ± . × − ( ± , ± , , ( ± , , ± , (0 , ± , ± . × − ( ± , ± , ± , ( ± , ± , ± , ( ± , ± , ± . × − ( ± , , , (0 , ± , , (0 , , ± . × − ( ± , , , (0 , ± , , (0 , , ± and finally, for N = 7, w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 0 . w = 5 . × − w = 5 . × − v = 0 , v = 1 , v = − , v = 2 , v = − v = 3 , v = − , v = 4 , v = − , v = 5 v = − , v = 6 , v = − c s = 1 . . (A.3)For the three dimensional case, up to fifth order, we can find the velocity vectors ~v i and weights w i by solving Eq. (26). The values are shown in Table 3. References
1. ¨Ottinger, H.C.,
Beyond Equilibrium Thermodynamics , Wiley, 2005.2. Zhenhua Chai and Baochang Shi, Applied Mathematical Modelling , 2050 (2008).3. Jackson John David, Electrodin´amica cl´asica , Editorial Alhambra S.A., first edition,1966. ovember 8, 2018 8:34 WSPC/INSTRUCTION FILELBPoissonApr25˙IJMPC
4. M. Ottaviani and F. Romanelli and R. Benzi and M. Briscolini and P. Santangelo andS. Succi, Physics of Fluids B: Plasma Physics , 67 (1990).5. Tomasi, J and Persico, M, Chemical Reviews , 2027 (1994).6. Honig, B and Nicholls, A, Science , 1144 (1995).7. Baker, NA and Sept, D and Joseph, S and Holst, MJ and McCammon, JA, Proc. ofthe Nat. Acad. of Sci. , 10037 (2001).8. Marques, M.A.L. and Gross, E.K.U., Annual Review of Physical Chemistry , 427(2004), PMID: 15117259.9. Sheng Meng and Efthimios Kaxiras, The Journal of Chemical Physics , 054110(2008).10. Gross, E.K.U. and Dreizler, R.M., Density Functional Theory , NATO ASI Series:Physics, Springer, 1995.11. R. Benzi and S. Succi and Vergassola, Phys. Rep. , 145 (1992).12. S. Chen and G. Doolen, Annu. Rev. Fluid Mech. , 329 (1998).13. Chen, Shiyi and Doolen, Gary D., Annual Review of Fluid Mechanics , 329 (1998).14. Ponce Dawson, S. and Chen, S. and Doolen, G. D., The Journal of Chemical Physics , 1514 (1993).15. Mendoza, M. and Mu˜noz, J. D., Phys. Rev. E , 056708 (2010).16. S. Succi and R. Benzi, Physica D
69 3-4 , 327 (1993).17. Mendoza, M. and Boghosian, B. M. and Herrmann, H. J. and Succi, S., Phys. Rev.Lett. , 014502 (2010).18. Mendoza, M. and Boghosian, B. M. and Herrmann, H. J. and Succi, S., Phys. Rev. D , 105008 (2010).19. Guo, Zhaoli and Zhao, T. S. and Shi, Yong, The Journal of Chemical Physics ,(2005).20. Jinku Wang and Moran Wang and Zhixin Li, Journal of Colloid and Interface Science , 729 (2006).21. Miki Hirabayashi and Yu Chen and Hirotada Ohashi, JSME International JournalSeries B Fluids and Thermal Engineering , 45 (2001).22. Xiaoyi He and Ning Li, Computer Physics Communications , 158 (2000).23. Baochang Shi and Bin Deng and Rui Du and Xingwang Chen, Computers and Math-ematics with Applications , 1568 (2008).24. Chikatamarla, Shyam S. and Karlin, Iliya V., Phys. Rev. E , 046701 (2009).25. Mendoza, M. and Karlin, I. and Succi, S. and Herrmann, H. J., Phys. Rev. D ,065027 (2013).26. Mendoza, M. and Succi, S. and Herrmann, H. J., Phys. Rev. Lett. , 096402 (2014).27. Chapman, S. and Cowling, T.G., The Mathematical Theory of Non-uniform Gases:An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion inGases , Cambridge Mathematical Library, Cambridge University Press, 1970.28. Bhatnagar, P. L. and Gross, E. P. and Krook, M., Phys. Rev. , 511 (1954).29. Grad, Harold, Communications on Pure and Applied Mathematics , 331 (1949).30. Sumesh P. Thampi and Santosh Ansumali and R. Adhikari and Sauro Succi, Journalof Computational Physics , 1 (2013).31. Rashmi Ramadugu and Sumesh P. Thampi and Ronojoy Adhikari and Sauro Succiand Santosh Ansumali, EPL (Europhysics Letters) , 50006 (2013).32. Hager, W.W., Applied numerical linear algebra , Prentice Hall PTR, 1988.33. Briggs, W.L. and Henson, V.E. and McCormick, S.F.,