sandbox/ghigo/src/test-particle/sphere-settling-large-domain-2.c

    Settling sphere in a large container at different Ga

    We investigate the settling of a sphere of diameter d in a large domain. This test case is governed by the Galileo number Ga = \frac{u d}{\nu}, where u = \sqrt{\lvert \frac{\rho_s}{\rho} - 1 \rvert d g} is the “steady” settling velocity.

    More precisely, we reproduce here 4 test cases taken from Uhlmann et al., 2004 (the article was redacted for an unknow reason), where the authors vary the Galileo number Ga to reproduce 4 different particle trajectory regimes:

    • \frac{\rho_{p}}{\rho}=1.5,\, Ga = 144 (case A): steady axisymmetric and vertical trajectory;
    • \frac{\rho_{p}}{\rho}=1.5,\, Ga = 178 (case B): steady oblique trajectory;
    • \frac{\rho_{p}}{\rho}=1.5,\, Ga = 190 (case C): low frequency oscillating oblique trajectory;
    • \frac{\rho_{p}}{\rho}=1.5,\, Ga = 250 (case D): chaotic or intermittent oblique trajectory.

    We solve here the 3D Navier-Stokes equations and describe the sphere using an embedded boundary.

    #include "grid/octree.h"
#include "../myembed.h"
#include "../mycentered.h"
#include "../myembed-particle.h"
#include "../myperfs.h"
#include "lambda2.h"
#include "view.h"

    Reference solution

    #define d       (1.) // Sphere diameter
#define grav    (2.) // Gravity
#define density (1.5)

#if GA // 
#define Ga   ((double) (GA))
#else // Ga = 178
#define Ga   (178) // Particle Reynolds number Re = ud/nu
#endif // GA

#define uref (sqrt (fabs ((density) - 1.)*(d)*(grav))) // Characteristic speed
#define tref ((d)/(uref)) // Characteristic time

    We also define the shape of the domain.

    #define sphere(x,y,z) (sq ((x)) + sq ((y)) + sq ((z)) - sq ((d)/2.))

void p_shape (scalar c, face vector f, coord p)
{
  vertex scalar phi[];
  foreach_vertex()
    phi[] = (sphere ((x - p.x),
		     (y - p.y),
		     (z - p.z)));
  boundary ({phi});
  fractions (phi, c, f);
  fractions_cleanup (c, f,
		     smin = 1.e-14, cmin = 1.e-14);
}

    We finally define the particle parameters.

    const double p_r = ((density)); // Ratio of solid and fluid density
const coord  p_i = {(p_moment_inertia_sphere ((d), (density))),
		    (p_moment_inertia_sphere ((d), (density))),
		    (p_moment_inertia_sphere ((d), (density)))}; // Particle moment of interia
const double p_v = (p_volume_sphere ((d))); // Particle volume
const coord  p_g = {0., -(grav), 0.};

    Setup

    We need a field for viscosity so that the embedded boundary metric can be taken into account.

    face vector muv[];

    We define the mesh adaptation parameters and vary the maximum level of refinement.

    #define lmin (9) // Min mesh refinement level (l=9 is 2pt/d)
#if LMAX // 11, 12, 13, 14
#define lmax ((int) (LMAX))
#else // 12
#define lmax (12) // Max mesh refinement level (l=12 is 16pt/d)
#endif // LMAX
#if CMAX
#define cmax ((((double) (CMAX))*1.e-3)*(uref))
#else // 1.e-2
#define cmax (1.e-2*(uref)) // Absolute refinement criteria for the velocity field
#endif // CMAX

    We also vary the CFL number.

    #if MYCFL
#define mycfl (((double) MYCFL)/100.)
#else
#define mycfl (0.5)
#endif // MYCFL

int main ()
{

    The domain is 256^3. It needs to be sufficiently big to allow for long settling times without the particle reaching the bottom.

      L0 = 256.;
  size (L0);
  origin (-L0/2., 0., -L0/2.);

    We set the maximum timestep.

      DT = 1.e-2*(tref);

    We set the tolerance of the Poisson solver.

      TOLERANCE    = 1.e-4;
  TOLERANCE_MU = 1.e-4*(uref);

    We initialize the grid.

      N = 1 << (lmin);
  init_grid (N);
  
  run();
}

    Boundary conditions

    We use no-slip boundary conditions.

    u.n[left] = dirichlet (0);
u.t[left] = dirichlet (0);
u.r[left] = dirichlet (0);

u.n[right] = dirichlet (0);
u.t[right] = dirichlet (0);
u.r[right] = dirichlet (0);

u.n[bottom] = dirichlet (0);
u.t[bottom] = dirichlet (0);
u.r[bottom] = dirichlet (0);

u.n[top] = dirichlet (0);
u.t[top] = dirichlet (0);
u.r[top] = dirichlet (0);

u.n[back] = dirichlet (0);
u.t[back] = dirichlet (0);
u.r[back] = dirichlet (0);

u.n[front] = dirichlet (0);
u.t[front] = dirichlet (0);
u.r[front] = dirichlet (0);

    We give boundary conditions for the face velocity to “potentially” improve the convergence of the multigrid Poisson solver.

    uf.n[left]   = 0;
uf.n[right]  = 0;
uf.n[bottom] = 0;
uf.n[top]    = 0;
uf.n[back]   = 0;
uf.n[front]  = 0;

    Properties

    event properties (i++)
{
  foreach_face()
    muv.x[] = (uref)*(d)/(Ga)*fm.x[];
  boundary ((scalar *) {muv});
}

    CFL condition

    event stability (i++)
{
  if (CFL > (mycfl))
    CFL = (mycfl);
}

    Initial conditions

    event init (i = 0)
{

    We set the viscosity field in the event properties.

      mu = muv;

    We use “third-order” face flux interpolation.

    #if ORDER2
  for (scalar s in {u, p, pf})
    s.third = false;
#else
  for (scalar s in {u, p, pf})
    s.third = true;
#endif // ORDER2

    We use a slope-limiter to reduce the errors made in small-cells.

    #if SLOPELIMITER
  for (scalar s in {u, p, pf}) {
    s.gradient = minmod2;
  }
#endif // SLOPELIMITER  
  
#if TREE

    When using TREE and in the presence of embedded boundaries, we should also define the gradient of u at the full center of cut-cells.

    #endif // TREE

    As we are computing an equilibrium solution when the particle reaches its settling velocity, we remove the Neumann pressure boundary condition which is responsible for instabilities.

    #if PNEUMANN // Non-zero Neumann bc for pressure
  for (scalar s in {p, pf}) {
    s.neumann_zero = false;
  }
#else // Zero Neumann bc for pressure
  for (scalar s in {p, pf}) {
    s.neumann_zero = true;
  }
#endif // PNEUMANN

    If the simulation is not restarted, we define the initial mesh and the initial velocity.

      if (!restore (file = "restart")) { // No restart

    We initialize the embedded boundary.

    We first define the particle’s initial position.

        p_p.y = (L0) - 5.*(d);
  
#if TREE

    When using TREE, we refine the mesh around the embedded boundary.

        astats ss;
    int ic = 0;
    do {
      ic++;
      p_shape (cs, fs, p_p);
      ss = adapt_wavelet ({cs}, (double[]) {1.e-30},
			  maxlevel = (lmax), minlevel = (1));
    } while ((ss.nf || ss.nc) && ic < 100);
#endif // TREE
  
    p_shape (cs, fs, p_p);
  }

    If the simulation is restarted, the proper restarting operations are performed in myembed-moving.h.

    }

    Embedded boundaries

    Adaptive mesh refinement

    #if TREE
event adapt (i++)
{
  adapt_wavelet ({cs,u}, (double[]) {1.e-2,(cmax),(cmax),(cmax)},
  		 maxlevel = (lmax), minlevel = (1));

    We do not need here to reset the embedded fractions to avoid interpolation errors on the geometry as the is already done when moving the embedded boundaries. It might be necessary to do this however if surface forces are computed around the embedded boundaries.

    }
#endif // TREE

    Restarts and dumps

    Every few characteristic time, we also dump the fluid data for post-processing and restarting purposes.

    #if DUMP
event dump_data (t += 5.*(tref))
{
  // Dump fluid
  char name [80];
  sprintf (name, "dump-level-%d-t-%g", (lmax), t/(tref));
  dump (name);

  // Dump particle
  char p_name [80];
  sprintf (p_name, "p_dump-level-%d-t-%g", (lmax), t/(tref));
  particle pp = {p_p, p_u, p_w, p_au, p_aw};
  struct p_Dump pp_Dump = {p_name, &pp};
  p_dump (pp_Dump);
}
#endif // DUMP

    Profiling

    #if TRACE > 1
event profiling (i += 20) {
  static FILE * fp = fopen ("profiling", "a"); // In case of restart
  trace_print (fp, 1); // Display functions taking more than 1% of runtime.
}
#endif // TRACE

    Outputs

    event coeffs (i++; t < 200.*(tref))
{
  char name1[80];
  sprintf (name1, "level-%d.dat", lmax);
  static FILE * fp = fopen (name1, "a"); // In case of restart
  fprintf (fp, "%d %g %g %d %d %d %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g\n",
	   i, t/(tref), dt/(tref),
	   mgp.i, mgp.nrelax, mgp.minlevel,
	   mgu.i, mgu.nrelax, mgu.minlevel,
	   mgp.resb, mgp.resa,
	   mgu.resb, mgu.resa,
	   p_p.x/(d), p_p.y/(d), p_p.z/(d),
	   p_u.x/(uref), p_u.y/(uref), p_u.z/(uref),
	   (Ga), fabs(p_u.y)*(d)/((uref)*(d)/(Ga))
	   );
  fflush (fp);

  double cell_wall = fabs (p_p.y - (d)/2.)/((L0)/(1 << (lmax)));
  if (cell_wall <= 1.)
    return 1; // stop
}

event logfile (t = end)
{
  fprintf (stderr, "%d %g %g %d %d %d %d %d %d %g %g %g %g %g %g %g %g %g %g %g %g\n",
	   i, t/(tref), dt/(tref),
	   mgp.i, mgp.nrelax, mgp.minlevel,
	   mgu.i, mgu.nrelax, mgu.minlevel,
	   mgp.resb, mgp.resa,
	   mgu.resb, mgu.resa,
	   p_p.x/(d), p_p.y/(d), p_p.z/(d),
	   p_u.x/(uref), p_u.y/(uref), p_u.z/(uref),
	   (Ga), fabs(p_u.y)*(d)/((uref)*(d)/(Ga))
	   );
  fflush (stderr);
}

    Snapshots

    event snapshots (t += 5.*(tref))
{
  scalar l2[];
  lambda2 (u, l2);
  boundary ({l2});

  char name2[80];

  clear();
  view (fov = 7, theta = 0, relative = false,
	tx = -(p_p.x)/(L0), ty = -(p_p.y + 25.*(d))/(L0),
	bg = {1,1,1},
	width = 400, height = 800);

  draw_vof ("cs", "fs");
  isosurface ("l2", -0.5, fc = {0.3,0.4,0.6}, lc = {0.3,0.4,0.6});
  sprintf (name2, "l2-5em1-level-%d-t-%g.png", lmax, t/(tref));
  save (name2);
  
  draw_vof ("cs", "fs");
  isosurface ("l2", -0.05, fc = {0.3,0.4,0.6}, lc = {0.3,0.4,0.6});
  sprintf (name2, "l2-5em2-level-%d-t-%g.png", lmax, t/(tref));
  save (name2);

  draw_vof ("cs", "fs");
  isosurface ("l2", -0.005, fc = {0.3,0.4,0.6}, lc = {0.3,0.4,0.6});
  sprintf (name2, "l2-5em3-level-%d-t-%g.png", lmax, t/(tref));
  save (name2);

  clear();
  view (fov = 2, theta = 0, relative = false,
	tx = -(p_p.x)/(L0), ty = -(p_p.y + 7.*(d))/(L0),
	bg = {1,1,1},
	width = 400, height = 800);

  draw_vof ("cs", "fs");
  isosurface ("l2", -0.5, fc = {0.3,0.4,0.6}, lc = {0.3,0.4,0.6});
  cells (n = {0, 0, 1}, alpha = (p_p.z));
  sprintf (name2, "l2-zoom-5em1-level-%d-t-%g.png", lmax, t/(tref));
  save (name2);
  
  draw_vof ("cs", "fs");
  isosurface ("l2", -0.05, fc = {0.3,0.4,0.6}, lc = {0.3,0.4,0.6});
  cells (n = {0, 0, 1}, alpha = (p_p.z));
  sprintf (name2, "l2-zoom-5em2-level-%d-t-%g.png", lmax, t/(tref));
  save (name2);

  draw_vof ("cs", "fs");
  isosurface ("l2", -0.005, fc = {0.3,0.4,0.6}, lc = {0.3,0.4,0.6});
  cells (n = {0, 0, 1}, alpha = (p_p.z));
  sprintf (name2, "l2-zoom-5em3-level-%d-t-%g.png", lmax, t/(tref));
  save (name2);
}

    Results on supercomputer (cedar) for CFL=0.5 and CMAX=1.e-2

    Settling velocity for Ga = 100

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.2:-1.0]
set title "Ga=100, CFL=0.5, cmax=1e-2"
plot "case-100/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-100/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-100/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-100/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-100/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-100/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-100/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-100/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 144

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.35:-1.0]
set title "Case A, Ga=144, CFL=0.5, cmax=1e-2"

# Steady settling values from Uhlamm et al, 2014
uAS = -1.292; # Case A, small domain
uAL = -1.285; # Case A, large domain
uAC15 = -1.2063; # Case AC-15, large domain
uAC36 = -1.2274; # Case AC-36, large domain

plot uAS w l lw 2 lc rgb "black" t "Ga=144 (AS), Uhlmann et al., 2014", \
     uAL w l lw 2 lc rgb "brown" t "Ga=144 (AL), Uhlmann et al., 2014", \
     uAC15 w l lw 2 lc rgb "purple" t "Ga=144 (AC-15), Uhlmann et al., 2014", \
     uAC36 w l lw 2 lc rgb "gray"   t "Ga=144 (AC-36), Uhlmann et al., 2014", \
     "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uAC15 - uAL)/uAL)*100. w l lw 2 lc rgb "purple" t "Ga=144 (AC-15), Uhlmann et al., 2014", \
     abs ((uAC36 - uAL)/uAL)*100. w l lw 2 lc rgb "gray"   t "Ga=144 (AC-36), Uhlmann et al., 2014", \
     "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-144/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 165

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.4:-1.0]
set title "Ga=165, CFL=0.5, cmax=1e-2"
plot "case-165/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-165/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-165/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-165/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-165/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-165/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-165/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-165/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 178

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.5:0]
set title "Case B, Ga=178, CFL=0.5, cmax=1e-2"

# Steady settling values from Uhlamm et al, 2014
uBS = -1.363; # Case B, small domain
uBL = -1.356; # Case B, large domain
uBC15 = -1.2434; # Case BC-15, large domain
uBC48 = -1.3067; # Case BC-48, large domain

plot uBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     uBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     uBC15 w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     uBC48 w l lw 2 lc rgb "gray" t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uBC15 - uBL)/uBL)*100. w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     abs ((uBC48 - uBL)/uBL)*100. w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:0.5]

# Steady horizontal values from Uhlamm et al, 2014
vBS = 0.1270; # Case B, small domain
vBL = 0.1245; # Case B, large domain
vBC15 = 0.2090; # Case BC-15, large domain
vBC48 = 0.1110; # Case BC-48, large domain

plot vBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     vBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     vBC15 w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     vBC48 w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "err(sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref})"
set yrange [0:200]
plot abs ((vBC15 - vBL)/vBL)*100. w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     abs ((vBC48 - vBL)/vBL)*100. w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the horinzontal velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]

# Steady angle values from Uhlamm et al, 2014
thetaBS = 5.323; # Case B, small domain
thetaBL = 5.225; # Case B, large domain

plot thetaBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     thetaBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.45:-1.25]
set yrange [*:*]
plot "case-178/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 190

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.5:0]
set title "Case C, Ga=190, CFL=0.5, cmax=1e-2"

# Average settling values from Uhlamm et al, 2014
uCS = -1.383; # Case C, small domain
uCL = -1.376; # Case C, large domain
uCC48 = -1.3233; # Case C-48 

plot uCS w l lw 2 lc rgb "black" t "Ga=190 (CS), Uhlmann et al., 2014", \
     uCL w l lw 2 lc rgb "brown" t "Ga=190 (CL), Uhlmann et al., 2014", \
     uCC48 w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uCC48 - uCL)/uCL)*100. w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:*]

# Average horizontal values from Uhlamm et al, 2014
vCS = 0.137; # Case C, small domain
vCL = 0.136; # Case C, large domain
vCC48 = 0.1201; # Case C-48 

plot vCS w l lw 2 lc rgb "black" t "Ga=190 (CS), Uhlmann et al., 2014", \
     vCL w l lw 2 lc rgb "brown" t "Ga=190 (CL), Uhlmann et al., 2014", \
     vCC48 w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "err(sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref})"
set yrange [0:200]
plot abs ((vCC48 - vCL)/vCL)*100. w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the horinzontal velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.45:-1.25]
set yrange [*:*]
plot "case-190/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 250

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.75:0]
set title "Case D, Ga=250, CFL=0.5, cmax=1e-2"

# Average settling values from Uhlamm et al, 2014
uDL = -1.4604; # Case D, large domain
uDC36 = -1.3946; # Case DC-36, large domain

plot uDL w l lw 2 lc rgb "brown" t "Ga=250 (DL), Uhlmann et al., 2014", \
     uDC36 w l lw 2 lc rgb "purple" t "Ga=250 (DC-36), Uhlmann et al., 2014", \
     "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uDC36 - uDL)/uDL)*100. w l lw 2 lc rgb "purple" t "Ga=250 (DC-36), Uhlmann et al., 2014", \
     "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:*]
plot "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.7:-1.3]
set yrange [*:*]
plot "case-250/CFL-0.50/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.50/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.50/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Results on supercomputer (cedar) for CFL=0.25 and CMAX=1.e-2

    Settling velocity for Ga = 100

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.2:-1.0]
set title "Ga=100, CFL=0.25, cmax=1e-2"
plot "case-100/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-100/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-100/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-100/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-100/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-100/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-100/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-100/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 144

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.35:-1.0]
set title "Case A, Ga=144, CFL=0.25, cmax=1e-2"

# Steady settling values from Uhlamm et al, 2014
uAS = -1.292; # Case A, small domain
uAL = -1.285; # Case A, large domain
uAC15 = -1.2063; # Case AC-15, large domain
uAC36 = -1.2274; # Case AC-36, large domain

plot uAS w l lw 2 lc rgb "black" t "Ga=144 (AS), Uhlmann et al., 2014", \
     uAL w l lw 2 lc rgb "brown" t "Ga=144 (AL), Uhlmann et al., 2014", \
     uAC15 w l lw 2 lc rgb "purple" t "Ga=144 (AC-15), Uhlmann et al., 2014", \
     uAC36 w l lw 2 lc rgb "gray"   t "Ga=144 (AC-36), Uhlmann et al., 2014", \
     "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uAC15 - uAL)/uAL)*100. w l lw 2 lc rgb "purple" t "Ga=144 (AC-15), Uhlmann et al., 2014", \
     abs ((uAC36 - uAL)/uAL)*100. w l lw 2 lc rgb "gray"   t "Ga=144 (AC-36), Uhlmann et al., 2014", \
     "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-144/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 165

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.4:-1.0]
set title "Ga=165, CFL=0.25, cmax=1e-2"
plot "case-165/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-165/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-165/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-165/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-165/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-165/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-165/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-165/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 178

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.5:0]
set title "Case B, Ga=178, CFL=0.25, cmax=1e-2"

# Steady settling values from Uhlamm et al, 2014
uBS = -1.363; # Case B, small domain
uBL = -1.356; # Case B, large domain
uBC15 = -1.2434; # Case BC-15, large domain
uBC48 = -1.3067; # Case BC-48, large domain

plot uBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     uBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     uBC15 w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     uBC48 w l lw 2 lc rgb "gray" t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uBC15 - uBL)/uBL)*100. w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     abs ((uBC48 - uBL)/uBL)*100. w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:0.5]

# Steady horizontal values from Uhlamm et al, 2014
vBS = 0.1270; # Case B, small domain
vBL = 0.1245; # Case B, large domain
vBC15 = 0.2090; # Case BC-15, large domain
vBC48 = 0.1110; # Case BC-48, large domain

plot vBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     vBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     vBC15 w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     vBC48 w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "err(sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref})"
set yrange [0:200]
plot abs ((vBC15 - vBL)/vBL)*100. w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     abs ((vBC48 - vBL)/vBL)*100. w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the horinzontal velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]

# Steady angle values from Uhlamm et al, 2014
thetaBS = 5.323; # Case B, small domain
thetaBL = 5.225; # Case B, large domain

plot thetaBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     thetaBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.45:-1.25]
set yrange [*:*]
plot "case-178/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 190

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.5:0]
set title "Case C, Ga=190, CFL=0.25, cmax=1e-2"

# Average settling values from Uhlamm et al, 2014
uCS = -1.383; # Case C, small domain
uCL = -1.376; # Case C, large domain
uCC48 = -1.3233; # Case C-48 

plot uCS w l lw 2 lc rgb "black" t "Ga=190 (CS), Uhlmann et al., 2014", \
     uCL w l lw 2 lc rgb "brown" t "Ga=190 (CL), Uhlmann et al., 2014", \
     uCC48 w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uCC48 - uCL)/uCL)*100. w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:*]

# Average horizontal values from Uhlamm et al, 2014
vCS = 0.137; # Case C, small domain
vCL = 0.136; # Case C, large domain
vCC48 = 0.1201; # Case C-48 

plot vCS w l lw 2 lc rgb "black" t "Ga=190 (CS), Uhlmann et al., 2014", \
     vCL w l lw 2 lc rgb "brown" t "Ga=190 (CL), Uhlmann et al., 2014", \
     vCC48 w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "err(sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref})"
set yrange [0:200]
plot abs ((vCC48 - vCL)/vCL)*100. w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the horinzontal velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.45:-1.25]
set yrange [*:*]
plot "case-190/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 250

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.75:0]
set title "Case D, Ga=250, CFL=0.25, cmax=1e-2"

# Average settling values from Uhlamm et al, 2014
uDL = -1.4604; # Case D, large domain
uDC36 = -1.3946; # Case DC-36, large domain

plot uDL w l lw 2 lc rgb "brown" t "Ga=250 (DL), Uhlmann et al., 2014", \
     uDC36 w l lw 2 lc rgb "purple" t "Ga=250 (DC-36), Uhlmann et al., 2014", \
     "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uDC36 - uDL)/uDL)*100. w l lw 2 lc rgb "purple" t "Ga=250 (DC-36), Uhlmann et al., 2014", \
     "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:*]
plot "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.7:-1.3]
set yrange [*:*]
plot "case-250/CFL-0.25/lmax-12/cmax-10e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-10e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-10e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Results on supercomputer (cedar) for CFL=0.25 and CMAX=1.e-3

    Settling velocity for Ga = 100

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.2:-1.0]
set title "Ga=100, CFL=0.25, cmax=1e-3"
plot "case-100/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-100/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-100/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-100/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-100/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-100/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-100/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-100/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-100/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-100/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 144

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.35:-1.0]
set title "Case A, Ga=144, CFL=0.25, cmax=1e-3"

# Steady settling values from Uhlamm et al, 2014
uAS = -1.292; # Case A, small domain
uAL = -1.285; # Case A, large domain
uAC15 = -1.2063; # Case AC-15, large domain
uAC36 = -1.2274; # Case AC-36, large domain

plot uAS w l lw 2 lc rgb "black" t "Ga=144 (AS), Uhlmann et al., 2014", \
     uAL w l lw 2 lc rgb "brown" t "Ga=144 (AL), Uhlmann et al., 2014", \
     uAC15 w l lw 2 lc rgb "purple" t "Ga=144 (AC-15), Uhlmann et al., 2014", \
     uAC36 w l lw 2 lc rgb "gray"   t "Ga=144 (AC-36), Uhlmann et al., 2014", \
     "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uAC15 - uAL)/uAL)*100. w l lw 2 lc rgb "purple" t "Ga=144 (AC-15), Uhlmann et al., 2014", \
     abs ((uAC36 - uAL)/uAL)*100. w l lw 2 lc rgb "gray"   t "Ga=144 (AC-36), Uhlmann et al., 2014", \
     "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(abs (($18 - uAL)/uAL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-144/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-144/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-144/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 165

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.4:-1.0]
set title "Ga=165, CFL=0.25, cmax=1e-3"
plot "case-165/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [*:*]
plot "case-165/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-165/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-165/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-165/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-165/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-165/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [*:-1]
set yrange [*:*]
plot "case-165/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-165/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-165/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 178

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.5:0]
set title "Case B, Ga=178, CFL=0.25, cmax=1e-3"

# Steady settling values from Uhlamm et al, 2014
uBS = -1.363; # Case B, small domain
uBL = -1.356; # Case B, large domain
uBC15 = -1.2434; # Case BC-15, large domain
uBC48 = -1.3067; # Case BC-48, large domain

plot uBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     uBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     uBC15 w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     uBC48 w l lw 2 lc rgb "gray" t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uBC15 - uBL)/uBL)*100. w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     abs ((uBC48 - uBL)/uBL)*100. w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(abs (($18 - uBL)/uBL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:0.5]

# Steady horizontal values from Uhlamm et al, 2014
vBS = 0.1270; # Case B, small domain
vBL = 0.1245; # Case B, large domain
vBC15 = 0.2090; # Case BC-15, large domain
vBC48 = 0.1110; # Case BC-48, large domain

plot vBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     vBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     vBC15 w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     vBC48 w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "err(sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref})"
set yrange [0:200]
plot abs ((vBC15 - vBL)/vBL)*100. w l lw 2 lc rgb "purple" t "Ga=178 (BC-15), Uhlmann et al., 2014", \
     abs ((vBC48 - vBL)/vBL)*100. w l lw 2 lc rgb "gray"   t "Ga=178 (BC-48), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vBL)/vBL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the horinzontal velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]

# Steady angle values from Uhlamm et al, 2014
thetaBS = 5.323; # Case B, small domain
thetaBL = 5.225; # Case B, large domain

plot thetaBS w l lw 2 lc rgb "black" t "Ga=178 (BS), Uhlmann et al., 2014", \
     thetaBL w l lw 2 lc rgb "brown" t "Ga=178 (BL), Uhlmann et al., 2014", \
     "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.45:-1.25]
set yrange [*:*]
plot "case-178/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-178/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-178/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 190

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.5:0]
set title "Case C, Ga=190, CFL=0.25, cmax=1e-3"

# Average settling values from Uhlamm et al, 2014
uCS = -1.383; # Case C, small domain
uCL = -1.376; # Case C, large domain
uCC48 = -1.3233; # Case C-48 

plot uCS w l lw 2 lc rgb "black" t "Ga=190 (CS), Uhlmann et al., 2014", \
     uCL w l lw 2 lc rgb "brown" t "Ga=190 (CL), Uhlmann et al., 2014", \
     uCC48 w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uCC48 - uCL)/uCL)*100. w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(abs (($18 - uCL)/uCL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:*]

# Average horizontal values from Uhlamm et al, 2014
vCS = 0.137; # Case C, small domain
vCL = 0.136; # Case C, large domain
vCC48 = 0.1201; # Case C-48 

plot vCS w l lw 2 lc rgb "black" t "Ga=190 (CS), Uhlmann et al., 2014", \
     vCL w l lw 2 lc rgb "brown" t "Ga=190 (CL), Uhlmann et al., 2014", \
     vCC48 w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "err(sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref})"
set yrange [0:200]
plot abs ((vCC48 - vCL)/vCL)*100. w l lw 2 lc rgb "purple" t "Ga=190 (CC-48), Uhlmann et al., 2014", \
     "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(abs ((sqrt($17*$17 + $19*$19) - vCL)/vCL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the horinzontal velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.45:-1.25]
set yrange [*:*]
plot "case-190/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-190/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-190/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    Settling velocity for Ga = 250

    reset
set terminal svg font ",14"
set key font ",14" top right spacing 1.1
set xlabel "t/t_{ref}"
set ylabel "u_{p,y}/u_{ref}"
set xrange [0:200]
set yrange [-1.75:0]
set title "Case D, Ga=250, CFL=0.25, cmax=1e-3"

# Average settling values from Uhlamm et al, 2014
uDL = -1.4604; # Case D, large domain
uDC36 = -1.3946; # Case DC-36, large domain

plot uDL w l lw 2 lc rgb "brown" t "Ga=250 (DL), Uhlmann et al., 2014", \
     uDC36 w l lw 2 lc rgb "purple" t "Ga=250 (DC-36), Uhlmann et al., 2014", \
     "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:18 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:18 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:18 w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the settling velocity $u_ (script){p,y}

    set ylabel "err(u_{p,y}/u_{ref})"
set yrange [0:10]
plot abs ((uDC36 - uDL)/uDL)*100. w l lw 2 lc rgb "purple" t "Ga=250 (DC-36), Uhlmann et al., 2014", \
     "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(abs (($18 - uDL)/uDL)*100) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Error for of the settling velocity $u_ (script){p,y}

    set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set yrange [0:*]
plot "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"

    Time evolution of the velocity $(script){u_{p,x}

    set ylabel "angle (deg)"
set yrange [0:*]
plot "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:(180./pi*atan (sqrt($17*$17 + $19*$19)/abs($18))) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Angle (in degrees) of the particle motion with respect to the vertical (script)

    Angle (in degrees) of the particle motion with respect to the vertical (script)

    set ylabel "dt/t_{ref}"
set yrange [1.e-6:*]
set logscale y
plot "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 2:3 w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 2:3 w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 2:3 w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the time step (script)

    Time evolution of the time step (script)

    set xlabel "x/d"
set ylabel "z/d"
set zlabel "y/d"
set xyplane 0
set xrange [-5:5]
set yrange [-5:5]
set zrange [0:*]
unset logscale y
splot "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "blue"      t "lmax=12", \
      "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "red"       t "lmax=13", \
      "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 14:16:(-($15-256)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Time evolution of the particle trajectory (script)

    Time evolution of the particle trajectory (script)

    set xtics -2,0.05,-1
set xlabel "u_{p,y}/u_{ref}"
set ylabel "sqrt(u_{p,x}^2 + u_{p,z}^2)/u_{ref}"
set xrange [-1.7:-1.3]
set yrange [*:*]
plot "case-250/CFL-0.25/lmax-12/cmax-1e-3/level-12.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "blue"      t "lmax=12", \
     "case-250/CFL-0.25/lmax-13/cmax-1e-3/level-13.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "red"       t "lmax=13", \
     "case-250/CFL-0.25/lmax-14/cmax-1e-3/level-14.dat" u 18:(sqrt($17*$17 + $19*$19)) w l lw 2 lc rgb "sea-green" t "lmax=14"
    Phase diagram of the vertical and horizontal velocities (script)

    Phase diagram of the vertical and horizontal velocities (script)

    References

    [Uhlmann2014]

    Markus Uhlmann and Jan Dušek. The motion of a single heavy sphere in ambient fluid: a benchmark for interface-resolved particulate flow simulations with significant relative velocities. International journal of multiphase flow, 59:221–243, 2014.