Entropy Balance and Dispersive Oscillations in Lattice Boltzmann Models
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Entropy Balance and Dispersive Oscillationsin Lattice Boltzmann Models
Dave Packwood ∗ Abstract
We conduct an investigation into the dispersive post-shock oscillations in the entropic lattice-Boltzmann method (ELBM). To this end we use a root finding algorithm to implement the ELBM whichdisplays fast cubic convergence and guaranties the proper sign of dissipation. The resulting simulationon the one-dimensional shock tube shows no benefit in terms of regularization from using the ELBM overthe standard LBGK method. We also conduct an experiment investigating of the LBGK method usingmedian filtering at a single point per time step. Here we observe that significant regularization can beachieved.
Keywords : Fluid dynamics, lattice Boltzmann, entropy balance, dispersive oscillations, numerical test,shock tube
AMS subject classifications.
The Lattice Boltzmann methods in their original form (see [2, 12]) do not guarantee the proper entropyproduction and may violate the Second Law. The proper entropy balance remains up to now a challengingproblem for many lattice Boltzmann models [14].The Entropic lattice Boltzmann method (ELBM) was invented first in 1998 as a tool for construction ofsingle relaxation time lattice Boltzman models which respect the H -theorem [9]. For this purpose, insteadof the mirror image with local equilibrium as reflection center, the entropic involution was proposed, whichpreserves the entropy value. Later, we call it the Karlin-Succi involution [7]. In 2000, it was reportedthat exact implementation of the Karlin-Succi involution (which keeps the entropy balance) significantlyregularizes the post-shock dispersive oscillations [1]. This regularization seems very surprising, because theentropic lattice BGK (ELBGK) model gives a second-order approximation to the Navier–Stokes equation(different proofs of that degree of approximation were given in [12] and [4]), and due to the Godunov theorem[6] linear second-order finite difference methods have to be non monotonic.Moreover, Lax [10] and Levermore with Liu [11], demonstrated that these dispersive oscillations areunavoidable in classical numerical methods. Schemes with precise control of entropy production, studiedby Tadmor with Zhong [13], also demonstrated post-shock oscillations. Of course, there remains some gapbetween methods with proven existence of dispersive oscillations, and ELBM. However, recently, the existenceof oscillations in the vicinity of the shock, at small values of viscosity for ELBM, was reported for Burgers’ ∗ Department of Mathematics, University of Leicester, ([email protected]). α ( α = 0): S ( f + α ( f ∗ − f )) = S ( f ) , (1.1)where S is entropy, f is a current distribution, and f ∗ is the corresponding equilibrium.In the first part of this paper we discuss a different numerical implementation of the ELBGK and conductan investigation into exactly what stabilization properties it exhibits.The choice of the method for solution of (1.1) should be very precise, and in Section 4 we describe acubically converging root finding algorithm. It is not sufficient to have high precision when we have averagedeviation of the current distribution f from the associated equilibrium f ∗ . For example, for solutions withshocks, it is usual for the distribution of this deviation to far from being exponential [5], and there appearpoints with deviation of several orders higher than the average. Moreover, it is sufficient to smooth a solutionat one point only. We demonstrate this in the second part of the paper. We select the lattice site with mostnonequilibrium f and regularize the field of nonequilibrium entropy at this point with 3-point median filter[5]. As a result, the dispersive oscillations drastically decrease. The Lattice Boltzmann method arises as a discretization of Boltzmann’s kinetic transport equation ∂f∂t + v · ∇ f = Q c ( f ) . (2.2)The population function f describes the distribution of single particles in the system and the collision integral Q c their interaction. Altogether (2.2) describes the behaviour of the system at the microscopic level. Byselecting a finite number of velocities v i , ( i = 1 ,...,n ) we create discrete approximation of the kinetic equationin velocity space. An appropriate choice of the velocities and time step discretizes space. For a time stepof δt = 1 the lattice can be created by unscaled space shifts of the velocities, and we get the fully discretelattice Botzmann gas: f i ( x + v i ,t + 1) = f i ( x,t ) + Q i (2.3)where the proper transition from continuous collision integral Q c ( f ) to its fully discrete form { Q i } is assumed.The simplest and the most common choice for the discrete collision integral Q i is the Bhatnagar-Gross-Krookoperator with over-relaxation Q i = αβ ( f ∗ i − f i ) . (2.4)For the standard LBGK method α = 2 and β ∈ [0 ,
1] (usually, β ∈ [1 / , β = 1 / local equilibrium f ∗ i and β = 1 (the mirror reflection ) returns the collision for a liquid at the zero viscosity limit. For a viscosity ν the parameter β is chosen by β = δt/ (2 ν + δt ). It should be noted that a collision integral such as (2.4) demands priorknowledge of a local equilibrium state for the given lattice.2 variation on the LBGK is the ELBGK [1]. In this case α is varied to ensure a constant entropycondition according to the discrete H -theorem. In general the entropy function is based upon the lattice andcannot always be found explicitly. However in the case of the simple one dimensional lattice with velocities v = ( − c, ,c ) and corresponding populations f = ( f − ,f ,f + ) an explicit Boltzmann style entropy function isknown [8]: S ( f ) = − f − log( f − ) − f log( f / − f + log( f + ) . (2.5)With knowledge of such a function α is found as the non-trivial root of the equation S ( f ) = S ( f + α ( f ∗ − f )) . (2.6)The trivial root α = 0 returns the entropy value of the original populations. ELBGK then finds the non-trivial α such that (2.6) holds. This version of BGK collision one calls entropic BGK (or EBGK) collision.Solution of (2.6) must be found at every time step and lattice site. Entropic equilibria (also derived fromthe H -theorem) are always used for ELBGK. H -theorem for LBMs In the continuous case the Boltzmann H the Maxwellian distribution maximizes entropy and therefore alsohas zero entropy production. In the context of lattice Boltzmann methods a discrete form of the H -theoremhas been suggested as a way to introduce thermodynamic control to the system [9].From this perspective the goal is to find an equilibrium state equivalent to the Maxwellian in the contin-uum which will similarly maximize entropy. Before the equilibrium can be found an appropriate H functionmust be known for a given lattice. These functions have been constructed in a lattice dependent fashion in[8], and H = − S with S from (2.5) is an example of a H function constructed in this way.Using equilbria derived from a H function with entropy considerations in mind leads to a thermodynam-ically correct LBM. This is easy to see in the case of the EBGK collision operator (2.4) with explicit localequilibrium. EBGK collision obvioulsly respect the Second Law (if β ≤ α that with β = 1 (inviscid fluid) would give zero entropy production, thereforemaking the position of zero entropy production the limit of any relaxation. For the fixed α used in theLBGK method it remains possible, particularly for low viscosity fluids, to relax past this point resulting innegative entropy production, violating the Second Law.Near to the zero-viscosity limit the LBGK method produces spurious oscillations around shockwaves.Apart from the thermodynamic benefits of using ELBGK it has been claimed [1] that ELBGK’s thermo-dynamic considerations act as a regularizer. This claim seems to be at odds with other numerical methodswhich respect the same thermodynamic laws as ELBGK. For example the results of Tadmoor and Zhong[13] for an entropy correct method display intensive post-shock oscillations. Furthermore it has been demon-strated [10, 11] that such dispersive oscillations are artifacts of the lattice rather than thermodynamic issues.As ELBGK clearly operates on exactly the same lattice as LBGK and other finite difference schemes itwarrants a deeper investigation into exactly how it achieves the regularization properties claimed. In order to investigate the stabilization properties of ELBGK it is necessary to craft a numerical methodcapable of finding the non-trivial root in (2.6). In this section we fix the population vectors f and f ∗ , and3re concerned only with this root finding algorithm. We recast (2.6) as a function of α only: F ( α ) = S ( f + α ( f ∗ − f )) − S ( f ) . (4.7)In this setting we attempt to find the non-trivial root r of (4.7) such that F ( r ) = 0. It should be noted thatas we search for r numerically we should always take care that the approximation we use is less than r itself.An upper approximation could result in negative entropy production.The following theorem gives cubic convergence order for a simple algorithm for finding the roots of aconcave function based on local quadratic approximations to the target function. Analogously to the casefor Newton iteration, the constant in the estimate is the ratio of third and first derivatives in the interval ofiteration. Theorem.
For a three times continously differentiable concave entropy function F ( α ) an iterative rootfinding method based on the zeros of a second order Taylor parabola has cubic convergence sufficiently closeto the root. Proof:
Assume that we are operating in a neighbourhood r ∈ N , in which F ′ is negative (as well of course F ′′ is negative). At each iteration the new estimate for r is the greater root of the parabola P , the secondorder Taylor polynomial at the current estimate, P ( α ) = F ( α n ) + ( α − α n ) F ′ ( α n ) + ( α − α n ) F ′′ ( α n )2 . (4.8)The Lagrange remainder form of the error is F ( α ) = F ( α n ) + ( α − α n ) F ′ ( α n ) + ( α − α n ) F ′′ ( α n )2 + ( α − α n ) F ′′′ ( γ n )6= P ( α ) + ( α − α n ) F ′′′ ( γ n )6 , where γ n lies between α n and α . Evaluating this at r we see that | P ( r ) | ≤ | α n − r | sup a ∈ N | F ′′′ ( a ) | . Now, using the mean value theorem, for some value b n ∈ [ α n +1 ,r ], | P ( r ) | = | P ( α n +1 ) + ( r − α n +1 ) P ′ ( b n ) | ≥ | ( r − α n +1 ) | inf b ∈ N | ( F ′ ( b )) | . Combining the last two equations we see that | ( r − α n +1 ) | ≤ C | α n − r | , where C = 16 sup a ∈ N | F ′′′ ( a ) | (cid:30) inf b ∈ N | F ′ ( b ) | . (cid:4) We use a Newton step to estimate the accuracy of the method at each iteration: | α n − r | ≈ | F ( α n ) /F ′ ( α n ) | . (4.9)In fact we use a convergence criteria based not solely on α but on α || f ∗ − f || , this has the intuitive appealthat in the case where the populations are close to the local equilibrium ∆ S = S ( f ∗ ) − S ( f ) will be small and4 very precise estimate of α is unnecessary. We have some freedom in the choice of the norm used and weselect between the standard L norm and the entropic norm. The entropic norm is defined as || f ∗ − f || f ∗ = − (( f ∗ − f ) , D S (cid:12)(cid:12) f ∗ ( f ∗ − f )) , where D S (cid:12)(cid:12) f ∗ is the second differential of entropy at point f ∗ , and ( x,y ) is the standard scalar product.The final root finding algorithm then is beginning with the LBGK estimate x = 2 to iterate using theroots of successive parabolas. If this first initiation step produces non-positive population, then the positivityrule [4] could be used (instead of the mirror image we choose the closest value of α which gives non-negativevalue of populations). The same regularization rule might be suggested if there exists no root we are lookingfor. In the tests described below, this situation never arose.We stop the method at the point, | α n − r | · || f ∗ − f || < ǫ. (4.10)To ensure that we use an estimate that is less than the root, at the point where the method has convergedwe check the sign of F ( α n ). If F ( α n ) > F ( α n ) < α n = α n − F ( α n ) F ′ ( α n ) . At each time step before we begin root finding we eliminate all sites with ∆
S < − . For these siteswe make a simple LBGK step. At such sites we find that round off error in the calculation of F by solutionof equation (1.1) can result in the root of the parabola becoming imaginary. We note that in such cases amirror image given by LBGK is effectively indistinct from the exact ELBGK collision.We now experimentally study the convergence of the method. The convergence of the bisection methodis presented for control. For the bisection method we calculate an initial estimate using the root of theparabola (4.8) with α = 2. Whichever side of the root this estimate is on, an estimate for the oppositeside can be found using a double length Newton step. We then have both an upper and lower estimate forthe root as required for the beginning of the bisection method. For this test ǫ is set to 10 − . . This is themaximum accuracy following the bound on ∆ S of 10 − due to the quadratic nature of F .For shock tube test using 800 lattice sites at the 400th iteration step (see detailed description in Section 5),Fig. 1 shows that the parabola based method required two iterations, but not more than two, at some pointsin a vicinity of shock. In other areas one iteration is sufficient for the desired accuracy. Across the wholelattice the entropic norm stipulates a slightly greater number of iterations in both methods. A standard experiment for the testing of LBMs is the one-dimensional shock tube problem. The latticevelocities used are v = ( − , , ρ ( x ) = (cid:26) , ≤ x ≤ , . , ≤ x ≤ . Initially all velocities are set to zero. We compare the ELBGK equipped with the parabola based root findingalgorithm using the entropic norm with the standard LBGK method using both standard polynomial and5
00 200 300 400 500 600 700 800−0.500.511.522.5 xn a
100 200 300 400 500 600 700 800−0.500.511.522.5 xn b
100 200 300 400 500 600 700 80002468 xn c
100 200 300 400 500 600 700 80002468 xn d
Figure 1: Iterations required for convergence to ǫ = 10 − . under (4.10) using ( a ) Parabola method with L norm; ( b ) Parabola method with entropic norm; ( c ) Bisection method with L norm; ( d ) Bisection methodwith entropic norm. 6
200 400 600 8000.40.60.81 x ρ a x ρ b x ρ c x ρ d x ρ e x ρ f Figure 2: Density profile of the simulation of the shock tube problem following 400 time steps using ( a )LBGK with polynomial equilibria [ ν = (1 / · − ]; ( b ) LBGK with entropic equilibria [ ν = (1 / · − ];( c ) ELBGK [ ν = (1 / · − ]; ( d ) LBGK with polynomial equilibria [ ν = 10 − ]; ( e ) LBGK with entropicequilibria [ ν = 10 − ]; ( f ) ELBGK [ ν = 10 − ].entropic equilibria. The polynomial equilibria are given in [2, 12]: f ∗− = ρ (cid:0) − u + 3 u (cid:1) , f ∗ = 2 ρ (cid:18) − u (cid:19) , f ∗ + = ρ (cid:0) u + 3 u (cid:1) . The entropic equilibria also used by the ELBGK are available explicitly as the maximum of the entropyfunction (2.5), f ∗− = ρ − u − p u ) , f ∗ = 2 ρ − p u ) , f ∗ + = ρ u − p u ) . Now following (2.3) the governing equations for the simulation are f − ( x,t + 1) = f − ( x + 1 ,t ) + αβ ( f ∗− ( x + 1 ,t ) − f − ( x + 1 ,t )) ,f ( x,t + 1) = f ( x,t ) + αβ ( f ∗ ( x,t ) − f ( x,t )) ,f + ( x,t + 1) = f + ( x − ,t ) + αβ ( f ∗ + ( x − ,t ) − f + ( x − ,t )) . From this experiment we observe no benefit in terms of regularization in using the ELBGK rather than thestandard LBGK method (Fig. 2). In both the medium and low viscosity regimes ELBGK fails to supressthe spurious oscillations found in the standard LBGK method.7
100 200 300 400 500 600 700 8000.30.40.50.60.70.80.911.1 x ρ a x ρ b Figure 3: Density profile of the simulation of the shock tube problem following 400 time steps using ( a )LBGK with entropic equilibria and one point median filtering [ ν = (1 / · − ]; ( b ) LBGK with entropicequilibria and one point median filtering [ ν = 10 − ].To explain previous results showing regularization by the ELBGK we note that in the collision integral(2.4) that α and β are composite. In this sense entropy production controlled by α and viscosity controlledby β are the same thing. A weak lower approximation to α would lead effectively to addition of entropyat the mostly far from equilibrium sites and therefore would locally increase viscosity. This numericalviscosity could, probably, explain the regularization and smoothing of the shock profile seen in some ELBGKsimulations. Finally we consider regularizing the LBGK method using median filtering at a single point. We follow theprescription detailed in [5]. First, at each time step, we locate the single lattice site x with the maximumvalue of ∆ S ( x ), and call this value ∆ S x . Secondly, we find the median value of ∆ S in the three nearestneighbours of x including itself, calling this value ∆ S med . Now instead of being updated using the standardBGK over-relaxation this single site is updated as follows: f − ( x,t + 1) = f ∗− ( x + 1 ,t ) + r ∆ S med ∆ S x ( f − ( x + 1 ,t ) − f ∗− ( x + 1 ,t )) ,f ( x,t + 1) = f ∗ ( x,t ) + r ∆ S med ∆ S x ( f ( x,t ) − f ∗ ( x,t )) ,f + ( x,t + 1) = f ∗ + ( x − ,t ) + r ∆ S med ∆ S x ( f + ( x − ,t ) − f ∗ + ( x − ,t )) . We observe that filtering a single point at each time step still results in a significant amount of regularization(Fig. 3).We also examine in each case the lattice site where the filtering is applied. The zero position is definedas the rightmost lattice site with ∆
S >
250 −200 −150 −100 −50 00102030405060708090100 xn a −250 −200 −150 −100 −50 00102030405060708090100 xn b
Figure 4: Distribution of median filtering sites relative to the position of the shock following 400 time stepsusing ( a ) LBGK with entropic equilibria and one point median filtering [ ν = (1 / · − ]; ( b ) LBGK withentropic equilibria and one point median filtering [ ν = 10 − ]. We present three main conclusions from this study.1. We do not find any evidence that maintaining proper balance of entropy regularize spurious oscillationsthe Lattice Boltzmann method. For ELBGK we confirm the conclusions of Lax [10] and Levermorewith Liu [11] that dispersive oscillations are unavoidable in numerical simulation of shocks.2. In order to clean up the parasite dispersive oscillations in the Lattice Boltzmann method it is nec-essary to filter the entropy in some way, so as to reduce the extremely-localised incidents of highnon-equilibrium entropy; see [5]. Previously reported smoothing of shocks must have been via theinadvertent introduction of numerical dissipation. (Perhaps, this conclusion could be extended to allknown regularisers of LBM, including those proposed by ourselves in [4].)3. For the 1D shock tube, one only needs to filter the entropy at one point per time step (usually very localto the shock), even at very low viscosity, in order to effectively eliminate the post-shock oscillation.We can expect that in 2D and 3D shocks it will be also necessary to filter nonequilibrium entropy insome local maxima points near the shock front only. The entropy filtering for non-entropic equilibriais possible [5] with use of the Kullback–Leibler distance from current distribution to equilibrium (therelative entropy).The Matlab code used to produce these results is provided in the appendix.
References [1] S. Ansumali, I. V. Karlin, Stabilization of the Lattice Boltzmann method by the H -theorem: A numericaltest, Phys. Rev. E, 62 (6):7999–8003, 2000.[2] R. Benzi, S. Succi, and M. Vergassola, The lattice Boltzmann-equation – theory and applications, Phys.Reports, 222:145–197, 1992. 93] B. M. Boghosian, P. J. Love, and J. Yepez, Entropic lattice Boltzmann model for Burgers equation,Phil. Trans. Roy. Soc. A, 362:1691–1702, 2004.[4] R. A. Brownlee, A. N. Gorban, and J. Levesley, Stability and stabilization of the lattice Boltzmannmethod, Phys. Rev. E, 75:036711, 2007.[5] R. A. Brownlee, A. N. Gorban, and J. Levesley, Nonequilibrium entropy limiters in lattice Boltzmannmethods, Physica A, 387 (2-3):385–406, 2008.[6] S. K. Godunov, A Difference Scheme for Numerical Solution of Discontinuous Solution of HydrodynamicEquations, Math. Sbornik, 47:271–306, 1959.[7] A. N. Gorban, Basic types of coarse-graining, In: A. N. Gorban, N. Kazantzis, I. G. Kevrekidis,H.-C. ¨Ottinger, and C. Theodoropoulos (eds.), Model Reduction and Coarse-Graining Approaches forMultiscale Phenomena, pages 117–176. Springer, Berlin-Heidelberg-New York, 2006. cond-mat/0602024.[8] I. V. Karlin, A. Ferrante, and H. C. ¨Ottinger, Perfect entropy functions of the lattice Boltzmann method,Europhys. Lett. 47:182–188, 1999.[9] I. V. Karlin, A. N. Gorban, S. Succi, and V. Boffi, Maximum entropy principle for lattice kineticequations, Phys. Rev. Lett., 81:6–9, 1998.[10] P. D. Lax, On dispersive difference schemes, Phys. D, 18:250–254, 1986.[11] C. D. Levermore and J.-G. Liu, Large oscillations arising in a dispersive numerical scheme, Physica D99:191–216, 1996.[12] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Oxford University Press, NewYork, 2001.[13] E. Tadmor and W. Zhong, Entropy stable approximations of Navier–Stokes equations with no artificialnumerical viscosity, J. Hyberbolic Differ. Equ., 3:529–559, 2006.[14] W.-A. Yong and L.-S. Luo, Nonexistence of H theorems for the athermal lattice Boltzmann modelswith polynomial equilibria, Phys. Rev. E, 67:051105.10 Matlab Code %LBM Lattice Boltzmann Method% LBM(MC,V,NX,TM,MOV,ESPIL,LIMIT,NORM) solves the shock tube problem% with number of lattice% sites NX with separation of one and viscosity V upto TM time steps of% length one using method% MC:% 0 − LBGK polynomial equilibria,% 1 − LBGK entropic equilibria,% 2 − ELBM Newton iterations,% 3 − ELBM Bisection method% and entropic limiter% LIMIT:% 0 − No limiter% 1 − Median filtering% Output of a movie of the simulation can be controlled with% MOV:% 0 − No movie,% 1 − With movie% The accuracy of the root finding for ELBM can be controlled with% EPSIL in all cases the root found will result in entropy production.% Convergence is partly based on | | feq − f | | , the choice of the norm is% controlled with% NORM:% 0 − L1 norm% 1 − Entropic Norm;function [rho, mov, convergence, relativeFiltering, alpha] = ...LBM(methodChoice,viscosity,nx,timeMax,MOV,ELBMEpsilon,Limiter,Norm)% Over − relaxation parameterbeta = 1/(2*viscosity+1);time = 1; %Start Timerho = LBMInitialise(nx); %Densities at each lattice site are intialised as%the standard shock tube problemu = 0; %Velocities at each lattice site are initialised as zeropopulations = LBMQuasiEquilibria(rho,u,methodChoice,nx); %Left moving,%central and right moving populations at each lattice site respectively are%initialised as the quasiequilibriarelativeFiltering = zeros(1,251);if (MOV == 1)mov = moviein(timeMax);x = 1:nx;elsemov = 0;endwhile time ≤ timeMaxpopulations = LBMPropagate(populations,nx); %Propagate the populations%keyboard rho, u] = LBMLatticeParameters(populations); %Density and velocity are%calculated using the populationspopequilibriums = LBMQuasiEquilibria(rho,u,methodChoice,nx); %New% quasiequilibria are found using new density and velocity values[alpha , convergence] = LBMEntropicParameter(populations, ...popequilibriums,methodChoice,nx,ELBMEpsilon,Norm);% Finding the non trivial root for constant entropy in the ELBM, for% normal LBGK alpha = 0[limiterSites,shockRelative] = LBMLimiterSites(populations, ...popequilibriums,nx,Limiter);relativeFiltering(250 − shockRelative) = ...relativeFiltering(250 − shockRelative) + 1;populations = LBMCollide(populations,popequilibriums,beta,alpha,nx, ...Limiter,limiterSites);% Populations are collided with the quasiequilibriaif( ¬ isreal(populations))keyboardreturnendif(MOV == 1) % If movie parameter is enabled record a frameunlimitedSites = setdiff(1:nx,limiterSites);plot(unlimitedSites,rho(unlimitedSites),'.',limiterSites, ...rho(limiterSites),'r*');axis([0 nx 0.3 1.3])mov(:,time) = getframe;endtime = time + 1; %Increment timeend%A function to intialise the lattice densities for the shock tube problemfunction rho = LBMInitialise(nx)rho = zeros(1,nx);half = floor(nx/2);rho(1:half) = 1;rho(half+1:nx) = 0.5;%A function to compute given a current density vector either the polynomial%or entropic quasiequilbria.function equilibria = LBMQuasiEquilibria(rho,u,methodChoice,nx)equilibria = zeros(3,nx);if ( methodChoice == 0) %Polynomial Equilibriaequilibria(1,:) = rho./6.*(1 − − lse %Entropic Equilibriaequilibria(1,:) = rho./6.*( − − − sqrt(1 + 3.*u.ˆ2));equilibria(3,:) = rho./6.*( 3.*u − −
1) = populations(1,2:nx);populations(3,2:nx) = populations(3,1:nx − − populations(1,:))./rho;% Function to find the constant entropy parameter for ELBM, in the case of% LBKG this is simply 2function [alpha, convergence] = LBMEntropicParameter(populations, ...popequilibriums,methodChoice,nx,ELBMEpsilon,Norm)%Choice of method:%0 − LBGK polynomial equilibria,%1 − LBGK entropic equilibria,%2 − ELBM Parabola iterations,%3 − ELBM Bisection methodif (methodChoice == 0)alpha = 2*ones(1,nx);convergence = zeros(1,nx);elseif (methodChoice == 1)alpha = 2*ones(1,nx);convergence = zeros(1,nx);elseif (methodChoice == 2)convergence = zeros(1,nx);alpha = 2*ones(1,nx);Sf = LBMEvaluateS(populations,popequilibriums,zeros(1,nx));SEquil = LBMEvaluateS(popequilibriums,popequilibriums,zeros(1,nx));remaining = find(abs(SEquil − Sf) > − lpha(remaining) = LBMInteriorParabola(populations(:,remaining),...popequilibriums(:,remaining),Sf(remaining),alpha(remaining));Spop = LBMEvaluateS(populations(:,remaining), ...popequilibriums(:,remaining),alpha(remaining));Sdiff = LBMEvaluateDiffS(populations(:,remaining), ...popequilibriums(:,remaining),alpha(remaining)); ∆ Alpha(remaining) = (Spop − Sf(remaining))./Sdiff;nowdone = find( popnorm(remaining).*abs( ∆ Alpha(remaining)) ... ≤ ELBMEpsilon);nowremaining = find( popnorm(remaining).* ...abs( ∆ Alpha(remaining)) > ELBMEpsilon);done = remaining(nowdone);remaining = setdiff(remaining,done);SDone = LBMEvaluateS(populations(:,done), ...popequilibriums(:,done),alpha(done)) − Sf(done);negFind = find(SDone < − ∆ Alpha(negativeEntropy);SDone = LBMEvaluateS(populations(:,done), ...popequilibriums(:,done),alpha(done)) − Sf(done);if(SDone < ¬ isreal(alpha))keyboardendif(iteration > − Sf) > − − Sf(remaining); eftSide = find(Spop > ≤ ∆ Alpha(remaining) = (Spop − Sf(remaining))./Sdiff;right(remaining(leftSide)) = alpha(remaining(leftSide)) ... − ∆ Alpha(remaining(leftSide));left(remaining(rightSide)) = alpha(remaining(rightSide)) ...+ 2* ∆ Alpha(remaining(rightSide));if(right < left)keyboardenditeration = 1;mid = left + (right − left)/2;while (isempty(remaining) == 0)mid(remaining) = left(remaining) ...+ (right(remaining) − left(remaining))/2;Sleft = S(populations(:,remaining), ...popequilibriums(:,remaining),left(remaining));Sright = LBMEvaluateS(populations(:,remaining), ...popequilibriums(:,remaining),right(remaining));Smid = LBMEvaluateS(populations(:,remaining), ...popequilibriums(:,remaining),mid(remaining));nowdone = find( popnorm(remaining).*(right(remaining) ... − left(remaining)) ≤ ELBMEpsilon);nowremaining = find( popnorm(remaining).*(right(remaining) ... − left(remaining)) > ELBMEpsilon);done = remaining(nowdone);remaining = setdiff(remaining,done);alpha(done) = left(done);convergence(done) = iteration;for itr = 1:length(remaining)if (Smid(nowremaining(itr)) > Sf(remaining(itr)))left(remaining(itr)) = mid(remaining(itr));elseright(remaining(itr)) = mid(remaining(itr));endenditeration = iteration + 1;if( ¬ isreal(alpha)) eyboardendif(iteration > − populations);S = − alphapop(1,:).*log(alphapop(1,:)) ... − alphapop(2,:).*log(alphapop(2,:)./4) ... − alphapop(3,:).*log(alphapop(3,:));% A function to implement the norms necessary to measure convergence of the% root findingfunction norms = LBMNorms(populations,popequilibriums,nx,nChoice);norms = zeros(1,nx);if (nChoice == 0)for j = 1:nxnorms(j) = norm(populations(:,j) − popequilibriums(:,j),1);endelseif(nChoice == 1)for j = 1:nxpoptemp = (populations(1,j) ... − popequilibriums(1,j)).ˆ2./popequilibriums(1,j);poptemp = poptemp + (populations(2,j) ... − popequilibriums(2,j)).ˆ2./popequilibriums(2,j);poptemp = poptemp + (populations(3,j) ... − popequilibriums(3,j)).ˆ2./popequilibriums(3,j);norms(j) = sqrt(poptemp);endend% A function to find an interior approximation to the root of the entropy% parabolafunction intPab = LBMInteriorParabola(populations,popequilibriums, ...STarget,alpha)SAlphaZero = LBMEvaluateS(populations,popequilibriums,alpha);SDashAlphaZero = LBMEvaluateDiffS(populations,popequilibriums,alpha);SDash2AlphaZero = LBMEvaluateDiff2S(populations,popequilibriums,alpha);intPab = [];for j = 1:length(alpha)intPab(j) = max(roots([0.5.*SDash2AlphaZero(j) SDashAlphaZero(j) ...SAlphaZero(j) − STarget(j)])) + alpha(j);end% A function to find the first derivative of the entropy unction DiffS = LBMEvaluateDiffS(populations,popequilibriums,alpha)alphapop = populations + (ones(3,1)*alpha).*(popequilibriums ... − populations);partDiff = popequilibriums − populations;DiffS = − partDiff(1,:).*(log(alphapop(1,:)) + 1) ... − partDiff(2,:).*(log(alphapop(2,:)./4) + 1) ... − partDiff(3,:).*(log(alphapop(3,:)) + 1);% A function to find the second derivative of the entropy.function Diff2S = LBMEvaluateDiff2S(populations,popequilibriums,alpha)alphapop = populations + (ones(3,1)*alpha).*(popequilibriums − populations);partDiff = popequilibriums − populations;Diff2S = − partDiff(1,:).ˆ2./(alphapop(1,:)) ... − partDiff(2,:).ˆ2./(alphapop(2,:)) ... − partDiff(3,:).ˆ2./(alphapop(3,:));% A function to detect an appropriate lattice site for limiting and to% measure it's position relative to the leading edge of the shock.function [Sites,shockRelative] = LBMLimiterSites(populations,...popequilibriums,nx,Limiter)shockRelative = 0;if (Limiter == 0)Sites = [];elseif (Limiter == 1)Sf = LBMEvaluateS(populations,popequilibriums,zeros(1,nx));Sfeq = LBMEvaluateS(popequilibriums,popequilibriums,zeros(1,nx));DeltaS = Sfeq − Sf;shockPos = max(find(DeltaS > − − Sites;end%Function for simple BGK collision of populations with quasiequilibria with% a limiter applied if necessaryfunction populations = LBMCollide(oldpopulations,popequilibriums,beta, ...alpha,nx,Limiter,limiterSites)unlimitedSites = setdiff(1:nx,limiterSites);populations(:,unlimitedSites) = oldpopulations(:,unlimitedSites) ...+ beta.*(ones(3,1)*alpha(unlimitedSites)) ....*( popequilibriums(:,unlimitedSites) ... − oldpopulations(:,unlimitedSites) );if (Limiter == 1)Sf = LBMEvaluateS(oldpopulations(:,(limiterSites − − − − ∆ S = Sfeq − Sf;Smed = median( ∆ S);coeff = sqrt(Smed/ ∆ S(2)); opulations(:,limiterSites) = popequilibriums(:,limiterSites) ...+ coeff.*(oldpopulations(:,limiterSites) ... − popequilibriums(:,limiterSites));endpopequilibriums(:,limiterSites));end