Venus and Earth are attracted to each other due to gravity. Image courtesy of universe today.

    Venus and Earth are attracted to each other due to gravity. Image courtesy of universe today.

    Two Liquid Planets in 2D

    On this page we look at the implementation of a body force due to gravity. The set-up closely mimics three aqua planets that float around in a two dimensional space that interact with each other via the virtue of gravity.


    To model the flow inside the liquid planets and the effect of gravity we use a continuum discription. Therefore we invoke the Navier-Stokes solver. To distinguish the planets from their vacuum surroundings, we use the “two-phase” scheme including surface tension.

    #include "navier-stokes/centered.h"
    #include "two-phase.h"
    #include "tension.h"
    #include "view.h"

    Because our planets are small in the vast emptiness of space, we need a resolution corresponding to 8192 grid cells in each direction. The number of planets (np) is set to three. Furthermore, we declare fields for a bodyforce (grav) that will be calculated from a gravity potential (G). The later requies obvious physically-consistent boundary conditions.

    int maxlevel = 13;
    int np = 3;
    face vector grav[];
    scalar G[];
    G[left]   = dirichlet (0);
    G[right]  = dirichlet (0);
    G[top]    = dirichlet (0);
    G[bottom] = dirichlet (0);

    The grid is initialized with size 1024\times 1024 a.u., we let the reader decide if that is short for astronomical units or arbitrary units. We also set the fluid properties of the two phases. To omit the most prominent computational challenges, a small offer is made with respect to the relevant dimensionless numbers that discribe the aforementioned water-vacuum system.

    int main() {
      N = 256;
      L0 = 1024;
      X0 = Y0 = -L0/2;
      rho1 = 1.;
      rho2 = 0.001;
      mu1 = 0.01;
      mu2 = 0.0;
      f.sigma = 0.05;
      a = grav;


    To mimic the details of the big bang, we initialize our planets randomly. Therefore, we set a random seed. Special care is required to make sure all processors are on the same page when using non-shared-memory parellization (e.g. MPI).

    event init (t = 0) {
      srand (time(NULL));
      double p[4]; // for xp, yp, R, U
      refine (level < 10 && fabs(x) < 20. && fabs(y) < 20.);
      for (int j = 1; j <= np; j++) {
    # if _MPI
        if (pid()==0) {
    # endif
          p[0] = 20.*noise();
          p[1] = 20.*noise();
          p[2] = 3.*fabs(noise()) + 1.5;
          p[3] = sign(noise())* (0.5 + 20./sq(p[2])*fabs(noise()));
    # if _MPI
        MPI_Bcast (&p, 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    # endif
        double xp = p[0];
        double yp = p[1];
        double R = p[2];
        double U = p[3];
        refine (level < maxlevel && sq(x - xp)+sq(y - yp) < sq(R + 1) &&
    	    sq(x - xp) + sq(y - yp)>sq(R - 1.));
        scalar f1[];
        fraction (f1, sq(R) - sq(x - xp) - sq(y - yp));
        foreach() {
          f[] += f1[];
          u.x[] += U*f1[] *y/(sqrt(sq(x) + sq(y)));
          u.y[] += U*f1[]*-x/(sqrt(sq(x) + sq(y)));;
        f[] = min(f[], 1.);
      boundary (all);

    The Gravity Body Force

    As mentioned before, a gravitational potential can be evaluated from the (mass) density field. The relation between the two is discribed by a Poisson equation that can conviniently be solved for with the dedicated solver. The body force is then simply, F_g \propto -\nabla G. Remarkably, when using a domain decomposition for parralel simulations, this method faciliates the enheritly non-local gravity force to be evaluated globally, well done Multigrid strategy!

    event acceleration (i++) {
      poisson (G, f); 
      boundary ({G});
        grav.x[] -= f[]* (G[] - G[-1])/Delta;

    Sanity test

    In order to check if the equations of motions are solved accurately and if the set-up is physically consistent, we evaluate the centre of mass of the system and write the result to a file. For a ‘realistic’ system this postion should only move along a straight path, since no external forces are applied. Again, i.e. for two planents in an infinite and frictionless domain.

    event centre_of_mass (t += 1.) {
      double xm = 0;
      double ym = 0;
      double w = 0;
      foreach (reduction(+:xm) reduction(+:ym) reduction(+:w)) {
        xm += x*f[]*sq(Delta);
        ym += y*f[]*sq(Delta);
        w += f[]*sq(Delta);
      static FILE *  fp = fopen ("pos.dat", "w");
      fprintf (fp, "%g\t%g\n", xm/w, ym/w);

    Grid Adaptation

    Rather unsophisticated, refinement is only based on the estimated discretization error in the representation of the planet-fraction field.

    event adapt (i++)
      adapt_wavelet({f}, (double []){0.01}, maxlevel, 6);


    We generate a movie.

    double min;
    event set_min (i = 0; i < 10; i+= 5.)
      min = statsf (G).min;
    event bviewer (t += 0.1; t < 50.) {
      view (fov = 3.);
      draw_vof ("f", filled = 1, fc = {0.9, 0.1, 0.9});
      squares ("G", linear = true, min = min , max = min/2.);
      save ("mov.mp4");


    The movie below shows the gravity potential and the liquid planets in magenta;


    Great succes? We check if the centre of mass has moved in a straight line or not.

     set xlabel 'x[a.u]'
     set ylabel 'y [a.u]'
     set key left top box 1
     set size ratio -1
     plot 'pos.dat' using 1:2 t 'Centre of mass'


    Appearently, the results are not as acurate as we would have liked them to be.