Cosmological dynamics

    Consider an otherwise empty Newtonian universe with small but heavy objects. The force of gravity seems important. For this purpose we use inertial-particles.h and the Poisson solver for the gravity potential \Phi.

    scalar Phi[];
    vector dPhi[];
    #include "poisson.h"
    mgstats mgP;

    The inertial particles definitions must be tuned for this special purpose. First, particles each have a mass m. Note that we may loose compatibility with tracer-particles.h is special care if not taken.

    #ifndef ADD_PART_MEM
    #define ADD_PART_MEM coord u; double m; long unsigned int tag;
    #include "inertial-particles.h"
    Particles cosmos;
    event defaults (t = 0) {
      foreach_dimension() {
        Phi[left]  = dirichlet (0);
        Phi[right] = dirichlet (0);

    The particles accelerate according to the gravity potential gradient.

    \displaystyle \mathbf{a}_p = -\nabla \Phi.

    mgstats get_Phi (void) {
      particle_boundary (cosmos);
      scalar density[];
        density[] = 0;
      foreach_particle_in(cosmos) {  // ;)
        Point point = locate (p().x,p().y,p().z);
        density[] += p().m;
        density[] /= dv();
      return poisson (Phi, density);
    coord p_acc (PA_inp inp) {
      particle p = inp.p;
      coord a;
      double xp = p.x; double yp = p.y; double zp = p.z;  
        a.x = -interpolate (dPhi.x, xp, yp, zp); 
      return a;
    event gravity (i++) {
      mgP = get_Phi();
      gradients ({Phi}, {dPhi});