Algorithm refinement for the stochastic Burgers' equation - Alejandro

18 Pages·2006·0.49 MB·English

ARTICLE IN PRESS
JournalofComputationalPhysicsxxx(2006)xxx–xxx
www.elsevier.com/locate/jcp
Algorithm refinement for the stochastic Burgers’ equation
John B. Bell, Jasmine Foo, Alejandro L. Garcia *
CenterforComputationalSciencesandEngineering,LawrenceBerkeleyNationalLaboratory,Berkeley,CA94720,USA
Received15December2005;receivedinrevisedform3August2006;accepted20September2006
Abstract
Inthispaper,wedevelopanalgorithmrefinement(AR)schemeforanexcludedrandomwalkmodelwhosemeanfield
behaviorisgivenbytheviscousBurgers’equation.ARhybridsusetheadaptivemeshrefinementframeworktomodela
systemusingamolecularalgorithmwheredesiredwhileallowingacomputationallyfastercontinuumrepresentationtobe
usedintheremainderofthedomain.Thefocusinthispaperistheroleoffluctuationsonthedynamics.Inparticular,we
demonstratethatitisnecessarytoincludeastochasticforcingterminBurgers’equationtoaccuratelycapturethecorrect
behaviorofthesystem.Theconclusionwedrawfromthisstudyisthatthefidelityofmultiscalemethodsthatcoupledis-
paratealgorithmsdependsontheconsistentmodelingoffluctuationsineachalgorithmandonacoupling,suchasalgo-
rithm refinement,that preservesthis consistency.
(cid:2)2006Elsevier Inc. Allrights reserved.
Keywords: Burgers’equation;Stochasticpartialdifferentialequations;Adaptivemeshrefinement;Algorithmrefinement;Hybridmethods;
Asymmetricexcludedrandomwalk
1. Introduction
Algorithm refinement (AR) is an emerging paradigm in the modeling and simulation of multiscale prob-
lems. Mathematical models use distinctly different representations for microscopic and macroscopic scales
with the corresponding algorithms echoing this disparity. Particle-based algorithms are a class of methods,
typically used to model the microscopic scale, that represent the physical system by discrete, interacting enti-
ties.These‘‘particles’’representanythingfromindividualatomstoparcelsoffluidtobacteriatoautomobiles.
Field-basedalgorithms,typicallyusedtomodelthemacroscopicscale,arederivedfrommodelsbasedprimar-
ily on partial differential equations with the physical system represented by continuum fields.
Algorithm refinement schemes (sometimes called ‘‘multi-algorithm hybrids’’) couple structurally different
computationalschemessuchasparticle-basedmolecularsimulationswithcontinuumpartialdifferentialequa-
tion (PDE) solvers.1 The general idea is to perform detailed calculations using an accurate but expensive
* Correspondingauthor.Tel.:[phone];fax:[phone].
E-mailaddress:[e-mail](A.L.Garcia).
1 NotethatothertypesofARhybridsexist(e.g.,couplingspectralanddiscretealgorithms[33]).
0021-9991/$-seefrontmatter (cid:2)2006ElsevierInc.Allrightsreserved.
doi:10.1016/j.jcp.2006.09.024
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys. (2006),doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
2 J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx
algorithm in a small region (or for a short time), and couple this computation to a simpler, less expensive
methodappliedtotherest.TheformulationofanARschemerequires:projectingfromthemicroscopicmodel
tomacroscopic;refiningfrommacroscopictomicroscopic;andhandshakingbetweenthetworepresentations
wheretheyarecoupled.Arelatedissueistheestablishmentof‘‘refinementcriteria’’thatspecifywhenamicro-
scopic representation is needed and when a macroscopic representation is sufficient. Examples of algorithm
refinement applied to fluid dynamics may be found in [15,17,29,34,35,38]; AR hybrids for interfacial propa-
gation are discussed in [27,28,30].
One aspect of multiscale modeling that has received insufficient attention is the presence of spontaneous
fluctuations at microscopic scales and their effect on the macroscopic scale. Accurate modeling of many
phenomena require the correct representation of the variances and correlations of fluctuations, specifically
when studying systems where the microscopic stochastics drive a macroscopic phenomenon. For physical
systems, the correct treatment of fluctuations is especially important for stochastic, nonlinear systems, such
as those undergoing phase transitions, nucleation, noise-driven instabilities and combustive ignition. In
these and related applications, the nonlinearities can exponentially amplify the influence of the
fluctuations.
Stochastic fluctuations in AR schemes have been investigated for two simple diffusive systems: linear
diffusion [3,37] and the quasi-linear train model [5]. For those parabolic problems, one finds that when
a particle algorithm is coupled to a deterministic continuum algorithm the variance of fluctuations is
reduced in the particle regime near the interface. The variance of fluctuations within the continuum regime
falls quickly away from the interface; however, variables, such as fluid velocity in the train model, that
have long-range correlations retain these correlations of fluctuations (though at reduced magnitude) within
the deterministic continuum region. Finally, stochastic continuum algorithms may be formulated such that
when coupled to particle schemes they correctly duplicate the physical fluctuations throughout the compu-
tational domain.
Our longer term goal is to extend the development of AR methods with fluctuations to an adaptive mesh
and algorithm method for the fluctuating compressible Navier–Stokes using a framework analogous to the
non-fluctuating CNS solver discussed in [17]. As a prelude to that extension, in the present work we develop
an AR method for Burgers’ equation that couples nonlinear hyperbolic waves and diffusion. For the particle
(microscopic)modelweconsidertheasymmetricexcludedrandomwalk(AERW)[25,32].Thehydrodynamic
(macroscopic)modelfortheAERWmodelistheviscousBurgers’equationwithstochasticforcing[8,9,12,20].
Inthefollowingsection,wedescribeindetailtheAERWmodel.Inthefollowingsection,weintroducethe
formofBurgers’equationthatrepresentsthehydrodynamiclimitoftheAERWmodelanddescribeadiscret-
ization ofthat equation basedon a second-order Godunov scheme. In Section4, we discussthe construction
ofthehybridmethodthatusesanoveralladaptivemeshrefinementframeworktodesignthecouplingbetween
microscopic and macroscopic models.Section 5 contains computational examples that illustrateand validate
the hybrid algorithm. As discussed in the concluding section, the numerical results demonstrate the impor-
tanceandchallengeofaccuratelymodelingfluctuationstosimulateandresolvebothmicroscopicandmacro-
scopic phenomena.
2. Asymmetric excluded random walk
2.1. Theory
The microscopic model for our system is an asymmetric excluded random walk. This model was selected
since it has been studied extensively by the statistical mechanics community [25,32], shown to be equivalent,
in the hydrodynamic limit, to the stochastic Burgers’ equation and found to exhibit a variety of interesting
phenomena (e.g., long-ranged correlations of non-equilibrium fluctuations). The AERW model is a system
of N random walker particles on a two-dimensional rectangular lattice of dimensions M ·M . Each site is
x y
denoted by a coordinate pair (x,y ) where j=1,...,M and k=1,...,M .
j k x y
Only one particle may occupy a site; the occupation number n(x,y)=1 (or =0) if a site is occupied (or
unoccupied). We choose the horizontal, or x-dimension of the lattice to correspond to the spatial domain
of the PDE and define the corresponding density,
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys.(2006), doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx 3
1 XMy
u ðxÞ¼ nðx;y Þ ð1Þ
p j M j k
y k¼1
so 06u (x)61. At equilibrium u is homogeneous and binomially distributed as the sum of M Bernoulli
p p y
random variables, each with probability U=N/M M of occupation. The mean and variance are
x y
Uð1(cid:2)UÞ
hu i¼U; hdu2i¼hu2i(cid:2)hu i2 ¼ : ð2Þ
p p p p M
y
In the local equilibrium approximation the equal-time correlation of fluctuations is [32]
1
hdu ðxÞdu ðxÞi¼ hu ðxÞið1(cid:2)hu ðxÞiÞd : ð3Þ
p i p j M p i p i i;j
y
Atanon-equilibriumsteadystatethevarianceismorecomplicatedduetolong-rangedcorrelationsoffluctu-
ations [31].
Particles on the lattice move between adjacent sites according to the asymmetric exclusion process. Each
particle waits a random time between moves with a mean free time of s. The next particle to move is drawn
at random by choosing a random site (x,y ); if the site is occupied then its particle is selected otherwise
j k
another random site is chosen.
The selected particle may move up, down, left or right to an adjacent site, according to the probabilities
assigned to the system. We take the particles to move horizontally or vertically with equal probability
(p =p =1/2) and the probabilities of moving up or down conditioned on vertical movement as equal
l M
(p =p =1/2).Asymmetryisintroducingbytakingunequalconditionalprobabilitiesforattemptingtomove
› fl
leftorright,thatis,p 6¼p withp +p =1.Oncetheparticleandmovedirectionarechosen,theparticle
! !
moves to the destination site, if unoccupied; if the destination site is occupied then the particle remains in
place. In either case the time is advanced and the entire process repeats.
2.2. Numerics
Givenaninitialdensitydistributionu(x),thelatticeisinitializedbyrandomlyfillingsites.Thedynamicsis
advancedbyrandomlychoosingparticlesandmovedirections,asdescribedabove.Inparticlesimulationsthe
physicaltimemaybeadvancedcontinuously(e.g.,event-driven dynamics)orintimeincrements(e.g.,molec-
ular dynamics) and either approach may be used for the AERW. For the former, the time between moves is
chosen asan exponentialrandom variablewithmean s/N,whereNisthenumberofparticles.For thelatter,
thenumberofmovesthatoccurduringatimeincrementDt isaPoissondistributedrandomvaluewithmean
p
l=NDt /s;ifDt (cid:3)s/NthentheprobabilityofamoveoccurringduringaparticletimestepDt isl+O(l2).
p p p
The lattice is periodic in y so that particles attempting to move up from row M move to the bottom
y
(first) row, provided it is unoccupied, with a similar definition for particles at the bottom row attempting
to move downward. If the x-direction is also periodic, then its treatment is analogous to the treatment
of periodicity in y.
Theother typeofboundaryconditionweconsider istheimpositionofDirichlet conditionsinx;inpartic-
ular, fixing particle densities, u and u at the left and right boundaries, respectively. These boundary condi-
L R
tions represent the occupation probabilities for each site on the boundary. We view the system as being
augmented with fictitious columns at j=0 and j=M +1 and with an effective total number of particles
x
N ¼N þu M þu M : ð4Þ
e L y R y
WethenviewtheAERWasoccurringontheenlargedlattice((M +2)·M )withprobabilistic‘‘virtual’’par-
x y
ticlesinthetwoboundarycolumns.WenotethatwithDirichletboundaryconditions,thenumberofparticles
in the system, N varies in time. Operationally these virtual particles enter the algorithm in two ways. First,
supposetheselectedparticlelocationforthenextmoveisintheleftboundarycolumn,say(x ,y );withprob-
0 k
abilityu thatsiteisconsideredoccupiedbyavirtualparticle.Iftheadjacentsite,(x ,y ),isunoccupiedthen
L 1 k
with probability p p a virtual particle moves to that destination, becoming a real particle. Similarly, if a
M !
particleattemptstojumpintotheleftboundaryfromaninteriorposition,thedestinationisunoccupiedwith
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys. (2006),doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
4 J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx
probability1(cid:2)u inwhichcasethejumpisacceptedandtheparticleremoved.Analogousrulesapplytothe
L
right boundary.
3. Burgers’ equation and continuum method
3.1. Theory
TheAERWmodeloftheprevioussectionisdefinedentirelyintermsofadiscretelattice.Inordertodefine
amacroscopicmodel,wespatiallyembedtheAERWmodelbyassigningaspatialwidthDx tothelatticesites
p
in the x direction. With this definition the hydrodynamic limit of the asymmetric excluded random walk
described in the previous section is the stochastic Burgers’ equation [20]:
o o
uðx;tÞ¼(cid:2) ffðuÞþdðuÞþgðx;tÞg; ð5Þ
ot ox
where g is a stochastic flux and
ou
fðuÞ¼cuð1(cid:2)uÞ; dðuÞ¼(cid:2)(cid:2) ð6Þ
ox
arethenonlinearadvectiveflux(soundspeed,c)andthediffusionflux(diffusionconstant,(cid:2)).Withthechange
of variable u0=(1(cid:2)2u)c this may be written in the more traditional form,
u0þu0u0 ¼(cid:2)u0 þ2cg ð7Þ
t x xx x
Note that variants of the stochastic Burgers’ equation, with different types of stochastic forcing, are com-
mon in the literature (e.g. [8]). Also note that there are other particle models, such as the Boghosian–Lev-
ermore cellular automaton, that also converge to a stochastic Burgers’ PDE in the hydrodynamic limit
[11,24].
The wave speed and diffusion constant are determined from the AERW parameters as
(cid:2) (cid:3)
2p Dx 1
c¼ $ p p (cid:2) ¼c ð2p (cid:2)1Þ ð8Þ
s ! 2 0 !
and
2p Dx2
(cid:2)¼ $ pp p ¼2c Dx p ð1(cid:2)p Þ; ð9Þ
s ! 0 p ! !
wherec =p Dx /sisthewavespeedforthecompletelyasymmetricwalk(p =0or1).Sincethewavespeed
0 M p !
f0(u) varies between +c and (cid:2)c on the range of u, we define a dimensionless cell Reynolds number as
jcDx j jp (cid:2)1j
Re ¼ p ¼ ! 2 ð10Þ
c (cid:2) p ð1(cid:2)p Þ
! !
thatcharacterizestherelativeimportanceofdiffusionandadvectionforthedynamicsatagivenmeshspacing.
Note that for p =1/2 the random walk is symmetric (pure diffusion) and Re =0; as p approaches 0 or 1
! c !
the random walk is unidirectional (pure advection) and Re goes to infinity.
c
The stochastic flux is a Gaussian white noise with zero mean and correlation
hgðx;tÞgðx0;t0Þi¼Aðx;tÞdðx(cid:2)x0Þdðt(cid:2)t0Þ; ð11Þ
wherethebracketsdenoteensembleaverage.Thenoiseamplitude,A(x,t),isrelatedthecorrelationofdensity
fluctuations; in the local equilibrium approximation,
hduðx;tÞduðx0;tÞi¼^uðx;tÞð1(cid:2)^uðx;tÞÞdðx(cid:2)x0Þ; ð12Þ
where ^u is the solution to the deterministic Burgers’ equation, that is,
^u ¼(cid:2)fð^uÞ þ(cid:2)^u : ð13Þ
t x xx
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys.(2006), doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx 5
From the above one finds,
Aðx;tÞ¼2(cid:2)^uð1(cid:2)^uÞ: ð14Þ
The noise amplitude may also be obtained from the continuum limit of the master equation for the AERW
[20].
3.2. Numerics
The stochastic Burgers’ equation may be simulated numerically by a variety of CFD algorithms, with the
choice guidedbytheapplication. Forexample, spectralmethods havebeendevelopedforhomogeneous, iso-
tropic turbulence (see [21,26] and references therein). For our AR hybrid we choose a cell-centered finite dif-
ference method, specifically a second-order Godunov scheme to calculate the hyperbolic flux and a simple
explicitpredictor-correctorcentereddifferenceschemetocomputethediffusionterm.See,forexample,Colella
[14]foradetaileddiscussionofsecond-orderGodunovmethods.Wedenotethespatialandtemporalgridsizes
asDxandDtanddenotebyuntheaverageofuincell-jattimen.Wenotethatfortheconstructionofahybrid
j
algorithm discussed later we will require that Dx be a multiple of Dx .
p
Thesecond-orderGodunovschemeconstructsalinearprofilewithineachcellwiththeslopesestimatedby
a higher-order finite difference approximation
(cid:2)un þ8un (cid:2)8un þun
u ¼ jþ2 jþ1 j(cid:2)1 j(cid:2)2: ð15Þ
x;j 12Dx
For advection-dominated problems a limiter istypically applied to these slopes; however,in the present con-
text, we are resolving atthe viscous length scale so no limiting isperformed. These slopes are used to predict
values at cell interfaces at the half-time level tn+1/2. In particular, we define
1
unþ1=2 ¼unþ ½Dx(cid:2)Dtmaxðf0ðunÞ;0Þ(cid:4)u ð16Þ
jþ1=2;‘ j 2 j x;j
and
1
unþ1=2 ¼un (cid:2) ½DxþDtminðf0ðunÞ;0Þ(cid:4)u ; ð17Þ
jþ1=2;r jþ1 2 j x;j
where f0=df/du. We then define the hyperbolic flux fnþ1=2 ¼fðunþ1=2Þ where unþ1=2 is the solution of the
jþ1=2 jþ1=2 jþ1=2
Riemann problem for u +f =0 along the ray x/t=0 with left and right states unþ1=2 and unþ1=2 ,
t x jþ1=2;‘ jþ1=2;r
respectively.
Thediffusionandstochasticfluxtermsareevaluatedusingapredictor-correctorscheme,treatingthehyper-
bolic flux terms as source terms. In particular, we first compute predicted values
p
ffiffiffi
Dt Dt 2Dt
up ¼un(cid:2) ðfnþ1=2(cid:2)fnþ1=2Þþ (cid:2)ðun (cid:2)2unþun Þþ ðgn (cid:2)gn Þ: ð18Þ
j j Dx jþ1=2 j(cid:2)1=2 ðDxÞ2 j(cid:2)1 j jþ1 Dx jþ1=2 j(cid:2)1=2
We then compute corrected values
up 1 Dt Dt pffi2ffiffiDt !
unþ1 ¼ j þ up(cid:2) ðfnþ1=2(cid:2)fnþ1=2Þþ (cid:2)ðup (cid:2)2upþup Þþ ðgp (cid:2)gp Þ : ð19Þ
j 2 2 j Dx jþ1=2 j(cid:2)1=2 ðDxÞ2 j(cid:2)1 j jþ1 Dx jþ1=2 j(cid:2)1=2
This can be rewritten as
Dt
unþ1 ¼un(cid:2) ðFn (cid:2)Fn Þ; ð20Þ
j j Dx jþ1=2 j(cid:2)1=2
where
1(cid:2) un (cid:2)un up (cid:2)up(cid:3) 1
Fn ¼fnþ1=2(cid:2) (cid:2) jþ1 j þ(cid:2) jþ1 j (cid:2)p ðgn þgp Þ ð21Þ
jþ1=2 jþ1=2 2 Dx Dx ffi2ffiffi jþ1=2 jþ1=2
is the total flux.
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys. (2006),doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
6 J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx
The stochastic flux for the asymmetric excluded random walk on an M ·M lattice is discretized as
x y
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
An;pþAn;p
gn;p ¼ j jþ1R: ð22Þ
jþ1=2 2DtM
y
where
An;p ¼2(cid:2)~un;pð1(cid:2)~un;pÞ; ~u¼minðmaxðu;0Þ;1Þ ð23Þ
j j j
and R is a Gaussian (normal) distributed random variable with zero mean and unit variance. Note that the
instantaneousfluctuatingvaluesareusedinplaceofthedeterministicvalue(seeEq.(14)),whichisaccurateas
long as the fluctuations remain small [18].
The scheme outlined aboveisstableprovided itsatisfiestimestep limitsforboththe hyperbolicand diffu-
sive terms. In particular, we require
Dtmaxjf0ðuÞj (cid:2)Dt 1
61 and 6 : ð24Þ
Dx Dx2 2
Wealsonotethatitispossibletouseasimplerversionoftheschemebasedonanexplicitfirst-ordertreatment
ofthediffusionterm.Forthemostpart,thesimplerversionprovidesreasonablepredictions;however,atlar-
gerDtthefirst-orderschemeover-predictsthevariationintheequilibriumsolutionbyapproximately5%,sug-
gesting that the temporal truncation error terms suppress the smoothing effect of the diffusion.
4. Algorithm refinement hybrid
In this section, we develop a hybrid algorithm refinement method that couples the AERW model intro-
duced in Section 2 with the stochastic Burgers’ equation algorithm in Section 3.
4.1. Basic construction
Philosophically,theconstructionofthehybridisbasedonthenotionthattheparticledescriptionprovides
amoreaccuraterepresentationofthesolutionthanthestochasticPDE.Thus,thebasicideaistorepresentthe
dynamicswiththecontinuummodelexceptinalocalizedregionwherehigher-fidelityparticlerepresentationis
required.
Our perspective in designing the algorithm follows the adaptive mesh and algorithm (AMAR) approach
introducedin[17].IncontrasttootherARapproaches(see,e.g.[15]),theAMARapproachmaintainsasolu-
tion of the macroscopic model over the entire domain (see Fig. 1). An error estimation criterion is used to
estimatewheretheimproved-representationoftheparticlemethodisrequired.Thatregion,whichcanchange
dynamically, is then ‘‘covered’’ with a particle patch. In this hierarchical representation, the solution is given
by the particle solution on the region covered by the particle patches and the continuum solution on the
remainder of the domain.
Thecouplingbetweentheparticleandcontinuumregionsusestheanalogofconstructsusedindeveloping
hierarchical adaptivemesh refinement algorithms.Forsimplicity, wewillassumethatthereisasinglerefined
patchandthatthemeshspacingforthecontinuumsolverisequaltothelatticespacingDx .Generalizationof
p
theapproachtoincludemultiplepatches(e.g.[39])andallowingthecontinuummeshtobeanintegralmulti-
ple of Dx (e.g. [4]) is fairly straightforward.
p
Integrationonthehierarchyisathreestepprocess.First,weintegratethecontinuumalgorithmfromtnto
tn+1,i.e.,foracontinuumstepDt.Theoldandnewstates,un andunþ1,areretaineduntiltheparticletimestep
j j
is complete. Continuum data at the edge of the particle patch is interpolated in time to provide Dirichlet
boundary conditions for the particle method.
Wehaveconsideredbothofthetime-evolutionschemesdiscussedinSection2.Fortheequaltimestepver-
sion,wechooseDt sothatDt =Dt/M foraspecifiedintegerM.WethenadvancetheparticlemethodbyM
p p t t t
steps until the particle and the continuum solver are at the same time. For the random time version of the
algorithm,theparticlemethodisadvancedbymoves,eachwitharandomtimeincrement,untilthenexttran-
sitionwouldadvancetheparticletimebeyondtn+1atwhichpointthetwosolutionsaresynchronized.Wenote
thatatthesynchronizationjuncturetheparticleandcontinuumsolutionsarenotquiteatthesametimelevel.
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys.(2006), doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx 7
ADVANCE
Particle Patch
Boundary Sites Boundary Sites
C
Continuum Grid B B
u
j
A
SYNCHRONIZE
D
E u E
j
Fig.1. SchematicofAERW/PDEhybrid.Advance:(A)advancecontinuumsolution;(B)setboundaryconditionsforAERWfromtime-
interpolatedPDEsolution;(C)advanceparticlesystembyAERW.Synchronize:(D)replaceoverlayingcontinuumvalueswithparticle
values;(E)resetPDEinterfacecellsbyrefluxing.
Forthemostpart,thishaslittleeffectonthecomputationalresults;however,itdoesleadtoerrorsofapprox-
imately 1% in the mean solution at equilibrium, which are not observed with the temporally synchronized
version.
4.2. Synchronization
The initial stage of the integration process essentially advances the macroscopic model separately with a
one-way coupling to the microscopic model by way of the Dirichlet boundary conditions. The macroscopic
model is not influenced by the microscopic model; the goal of the synchronization process is to correct the
macroscopic solution to reflect the effect of the microscopic model as though the integration were tightly
coupled.
Therearetwocomponentstothesynchronizationprocess.First,ontheregioncoveredbytheparticlerep-
resentation we replace the continuum solution obtained from the SPDE discretization by the more accurate
particle representation, i.e, set
1 XMy
unþ1 ¼ nðx;y Þ ð25Þ
j M j k
y k¼1
for each cell covered by the particle patch. Second, the continuum cells immediately adjacent to the particle
region,whichsuppliedboundarydatafortheparticleregionduringitsadvance,arecorrectedby‘‘refluxing’’.
Specifically,supposetheleftboundaryoftheparticlepatchoccursincellJ+1.ThevalueincontinuumcellJ
was updated with the continuum scheme
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys. (2006),doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
8 J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx
Dt
unþ1 ¼un (cid:2) ðFn (cid:2)Fn Þ ð26Þ
J J Dx Jþ1=2 J(cid:2)1=2
withthefluxesFncomputedfromthecontinuumvalues(seeEq.(20)).Themicroscopicallycorrectfluxisgiven
by the net number of particles moving across edge J+1/2 rather than by the continuum flux Fn . To per-
Jþ1=2
formtherefluxingcorrectionwemonitorthenumberofparticles,N! andN ,thatmoveintoandoutof
Jþ1=2 Jþ1=2
the particle region, respectively, across the continuum/particle interface at edge J+1/2. We then correct the
continuum solution as
Dt N! (cid:2)N
u0nþ1 ¼unþ1þ Fn (cid:2) Jþ1=2 Jþ1=2: ð27Þ
J J Dx Jþ1=2 DxM
y
This update effectively replaces the continuum flux component of the update to unþ1 on edge J+1/2 by the
J
fluxofparticlesthroughtheedge.Theuseofparticlefluxesisseamlessattheinterfacebecausethecontinuum
fluxisalreadystochastic(seeEqs.(21)and(22)).Ananalogousrefluxingstepoccursinthecelladjacenttothe
right-handboundaryoftheparticleregion.Finally,notethatthissynchronizationprocedureguaranteesexact
conservationofintegrateddensity.Thetechnicaldetailsofrefluxinginhigherdimensions(e.g.,thetreatment
of corners) are discussed in Garcia et al. [17].
4.3. Refinement criterion and regridding
TheAMARframeworkallowsustodynamicallychangethelocationoftheparticleregion.Therearesev-
eral possible strategies for designing refinement criteria. For the examples described in Section 5.3, we will
focus on criteria that identify cells where the solution has a large gradient characteristic of a viscous shock
profile. A straightforward measurement of the local gradient of the solution (e.g., (u (cid:2)u)/Dx) is not ade-
j+1 j
quate since the inherentfluctuations could trigger refinement even at equilibrium. What isneeded is a robust
measurethatidentifiesviscousshockswithoutgeneratingsubstantial‘‘falsepositives’’leadingtounnecessary
refinement. To this end we define a regional gradient using
" #
1 1 XS 1 XS
D ¼ u (cid:2) u ; ð28Þ
j SDx S jþi S j(cid:2)ði(cid:2)1Þ
i¼1 i¼1
where thestencil size Sisspecified; we takeS=4 inthecomputations in thefollowingsection.From Eq.(3)
onemayeasilyestimatetheexpectedstandarddeviationrofDresultingfromequilibriumfluctuationsandset
atoleranceofCr,whereCisaconstant.Toestimatewheretoplacetheparticleregionintheadaptivecodewe
computetheregionalgradientateachpoint;ifjDjexceedsthetolerancelevel,thencellsjandj+1aretagged.
j
Since we restrictourselves toa single particle patch, thelargestinterval containing all tagged cellsisthen the
newparticleregion.Ifmultiplepatcheswereallowed thentagged cellswouldbecollected toformparticlere-
gions;techniquesforcollecting tagged cellsinanoptimalmanner arewell established inthemesh refinement
literature [7].
Once the new particle region has been identified, it must be initialized. For continuum cells that were
already in the particle region, we simply retain the distribution of occupied sites. For cells that were not in
the previous particle region, we use the continuum density to compute N, the desired number of particles
f
forfillingacolumn.ThesimplewaytodothisistotakeN ¼unM randomlyroundedtothenearestinteger;
f j y
an alternative approach would be to fill randomly each site with probability un. We use the former approach
j
sinceitpreservesconservationoftotaldensity(towithinquantizationrounding).Wenotetheregriddingalgo-
rithmdoesnotneedtobedoneeverystep.SimpleestimatesbasedeitheronCFLconsiderationsorestimates
of discrete traveling wave velocities can be used to determine how often to regrid [6].
5. Computational examples
This section presents a series of computational examples, of progressively increasing sophistication, that
demonstrate the accuracy and effectiveness of the algorithm refinement hybrid. We consider four numerical
schemes:theasymmetricexcludedrandomwalk(AERW)fromSection2;thestochasticPDE(SPDE)forBur-
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys.(2006), doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx 9
gers’ equation from Section 3 and two algorithm refinement hybrids from Section 4. The first hybrid couples
theAERWandSPDE,withtheparticleschemeinasinglepatchwithinthesystem.Thesecondhybridissim-
ilar but without a stochastic flux, that is, using a deterministic PDE (DPDE). Both fixed-patch and adaptive
hybrids are considered as well as a handful of minor variants, discussed below.
Mean
0.52
0.515
0.51
0.505
0.5
0.495
0.49
0.485
0.48
–0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 0.2 0.3 0.4 0.5
x 10–3 Variance
1.9
1.8
1.7
1.6
1.5
1.4
–0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 0.2 0.3 0.4 0.5
x 10–5 Correlation
5
Hybrid w/ dt= 0.05
Hybrid w/ dt = 0.1
Hybrid w/ dt = 0.2
AERW
SPDE w/dt = 0.05
0
–5
–0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 0.2 0.3 0.4 0.5
x
Fig.2. MeanÆuæ,varianceÆdu2(x)æandcenterpointcorrelationÆdu(x)du(0)æversusxforasystematequilibrium(u =u =0.5).Lines
L R
are:SPDE/AERWhybridwithDt=0.05(reddasheddot);SPDE/AERWhybridwithDt=0.1(solidred);SPDE/AERWhybridwith
Dt=0.2(reddotted);AERW(green);SPDEwithDt=0.05(solidblue);predictedvarianceÆuæ(1(cid:2)Æuæ)/M (dashedblack).Notethatfor
y
thecorrelationÆdu(x)du(0)æthelargespikeatx=0isomittedfromtheplot.
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys. (2006),doi:10.1016/j.jcp.2006.09.024
ARTICLE IN PRESS
10 J.B.Belletal./JournalofComputationalPhysicsxxx(2006)xxx–xxx
The following parameter values are common to the simulations in all the examples: M =150,
y
Dx=Dx =0.01, s=1, p =1/2, p =p =1/2, and c =0.01. The particle time step is chosen to be small
p M › fl 0
(Dt (cid:5)s/(75N)) so the probability of a move occurring during a time step is taken as NDt /s. For a system
p p
length of L the viscous relaxation time is T =L2/(cid:2); in our simulations T =O(104). For the simulations of
(cid:2) (cid:2)
steadystatesthesystemisinitializednearthefinalstateandallowedtorelaxforatimethatislongcompared
to the relaxation time (typically for >100T ) before taking samples.
(cid:2)
5.1. Equilibrium state
First, we consider the simplest scenario, a system at the equilibrium state with equal, fixed density at
x=(cid:2)0.5and0.5.Theprobabilityofmovingtotherightisp =0.55,correspondingtoacellReynoldsnum-
!
ber of Re =0.20 (weakly hyperbolic). The single-algorithm simulations (AERW and Burger’s SPDE) have
c
100 sites or grid points in the x-direction. The AR hybrids introduce a fixed particle patch at the center of
the system between x=(cid:2)0.1 and x=0.1 with M =20. The hybrid simulation is performed at three contin-
x
uumtimestepsizes:Dt=0.05,0.1,and0.2.Sinceincrementaltimesteppingisusedintheparticlealgorithm,
to keep Dt =Dt/M (cid:3)Dt we take three corresponding values: M =8000, 14,000, and 28,000. Each simula-
p t t
tion is run to a final time of T=2·107, which corresponds to N =T/Dt continuum time steps (e.g.,
t
N =4·108 for the smallest Dt).
t
TypicalresultsfromthevariousnumericalschemesareshowninFig.2wherethemean,Æuæ;variance,Ædu2æ;
andcorrelation,Ædudu0æofdensityareplottedversusposition.Thesethreequantitiesareestimatedfromsam-
ples as
x 10–5 Correlation
1
SPDE: 1
SPDE: 2
AERW
0
–1
–2
–0.5 –0.4 –0.3 –0.2 –0.1 0 0.1 0.2 0.3 0.4 0.5
x
Fig. 3. Center point correlation versus x for a system at equilibrium (u =u =Æuæ=0.5) with the SPDE method using mean and
L R
instantaneous solution for the noise amplitude. Lines are SPDE: 1 using the mean (blue); SPDE: 2 using the instantaneous solution
(magenta);AERW(green).TimestepisDt=0.05.
Pleasecitethisarticleinpressas:J.B.Bellet al.,AlgorithmrefinementforthestochasticBurgers’equation, J.Comput.
Phys.(2006), doi:10.1016/j.jcp.2006.09.024