brownian.c

Brownian motion

Brownian motion

A hundred million particles take a random walk. We compare their behaviour as a group against the diffusion of a scalar field.

set xlabel 'x'
set ylabel 'pdf'
set grid 
plot 'log' w l lw 2 t 'Initial Particle PDF' ,\
     'log' u 1:3 t 'Scalar (t = 0)' ,\
     'out' w l lw 2 t 'Final particle PDF' ,\
     'out' u 1:3 t 'Scalar (t = t_{end})'

The statistics of a random walk appear as diffusion (script)

#include "grid/multigrid1D.h"
#include "particle.h"
#include "diffusion.h"
#include "run.h"

Particles Brownian;
int np = 1e8;

scalar a[];

double Diff = 0.25, tend = 5;

int main() {
  periodic (left);
  L0 = 12.;
  X0 = -L0/3.;
  N  = 128;
  run();
}

event init (t = 0) {
  DT = 0.05;
  Brownian = new_particles (np);

Initially, the particles are placed `normally’, following the Box–Muller transform.

  foreach_particle () {
    double a = rand()/(double)(RAND_MAX);
    double b = rand()/(double)(RAND_MAX);
    p().x = sqrt(-2*log(a))*cos(2*pi*b);
  }
  scalar s[];
  particle_pdf (Brownian, s);
  foreach() {
    a[] = exp(-sq(x)/2.)/sqrt(2*pi);
    fprintf (stderr, "%g %g %g\n", x, s[], a[]);
  }
}

event integrate (i++) {
  dt = dtnext (DT);
  
  const face vector mu[] = {Diff, Diff, Diff};
  diffusion (a, dt, mu);

A random walk is performed with Einstein’s step size.

  random_step (Brownian, sqrt(2.*Diff*dt));
}

event stop (t = tend) {
  particle_boundary (Brownian); //BCs are applied *after* the simulation
  scalar s[];
  particle_pdf (Brownian, s);
  foreach()   
    printf ("%g %g %g\n", x, s[], a[]);
  return 1;
}