/**
# Circular droplet in equilibrium

This is the classical "spurious" or "parasitic currents" test case
discussed in [Popinet, 2009](/src/references.bib#popinet2009).

We use the Navier--Stokes solver with VOF interface tracking and
surface tension. */

#define JACOBI 1

#include "navier-stokes/centered.h"
#include "vof.h"
#include "tension.h"

/**
The interface is represented by the volume fraction field *c*. */

scalar c[], * interfaces = {c};

/**
The diameter of the droplet is 0.8. The density is constant (equal to
unity by default), and the viscosity is defined through the Laplace
number
$$
La = \sigma\rho D/\mu^2
$$
with $\sigma$ set to one. The simulation time is set to the
characteristic viscous damping timescale. */

#define DIAMETER 0.8
#define MU sqrt(DIAMETER/LAPLACE)
#define TMAX (sq(DIAMETER)/MU)

/**
We will vary the number of levels of refinement (to study the
convergence), the Laplace number and *DC* a convergence parameter
which measures the variation in volume fraction between successive
timesteps (to evaluate whether we are close to a steady solution). */

int LEVEL;
double LAPLACE;
double DC = 0.;
FILE * fp = NULL;

int main() {
  
  /**
  We neglect the advection terms and vary the Laplace, for a constant
  resolution of 5 levels. */

  TOLERANCE = 1e-6 [*];
  stokes = true;
  c.sigma = 1;
  LEVEL = 5;
  N = 1 << LEVEL;
  for (LAPLACE = 120; LAPLACE <= 12000; LAPLACE *= 10)
    run();

  /**
  We now fix the Laplace number and look for stationary solutions
  (i.e. the volume fraction field is converged to within 1e-10) for
  varying spatial resolutions. */

  LAPLACE = 12000; DC = 1e-10;
  for (LEVEL = 3; LEVEL <= 7; LEVEL++) 
    if (LEVEL != 5) {
      N = 1 << LEVEL;
      run();
    }
}

/**
We allocate a field to store the previous volume fraction field (to
check for stationary solutions). */

scalar cn[];

event init (i = 0) {

  /**
  We set the constant viscosity field... */

  const face vector muc[] = {MU,MU};
  mu = muc;

  /**
  ... open a new file to store the evolution of the amplitude of
  spurious currents for the various LAPLACE, LEVEL combinations... */

  char name[80];
  sprintf (name, "La-%g-%d", LAPLACE, LEVEL);
  if (fp)
    fclose (fp);
  fp = fopen (name, "w");

  /**
  ... and initialise the shape of the interface and the initial volume
  fraction field. */
  
  fraction (c, sq(DIAMETER/2) - sq(x) - sq(y));
  foreach()
    cn[] = c[];
}

event logfile (i++; t <= TMAX)
{
  /**
  At every timestep, we check whether the volume fraction field has
  converged. */
  
  double dc = change (c, cn);
  if (i > 1 && dc < DC)
    return 1; /* stop */

  /**
  And we output the evolution of the maximum velocity. */

  scalar un[];
  foreach()
    un[] = norm(u);
  fprintf (fp, "%g %g %g\n",
	   MU*t/sq(DIAMETER), normf(un).max*sqrt(DIAMETER), dc);
}

event error (t = end) {
  
  /**
  At the end of the simulation, we compute the equivalent radius of
  the droplet. */

  double vol = statsf(c).sum;
  double radius = sqrt(4.*vol/pi);

  /**
  We recompute the reference solution. */
  
  scalar cref[];
  fraction (cref, sq(DIAMETER/2) - sq(x) - sq(y));
  
  /**
  And compute the maximum error on the curvature *ekmax*, the norm of
  the velocity *un* and the shape error *ec*. */
  
  double ekmax = 0.;
  scalar un[], ec[], kappa[];
  curvature (c, kappa);
  foreach() {
    un[] = norm(u);
    ec[] = c[] - cref[];
    if (kappa[] != nodata) {
      double ek = fabs (kappa[] - (/*AXI*/ + 1.)/radius);
      if (ek > ekmax)
	ekmax = ek;
    }
  }
  
  /**
  We output these on standard error (i.e. the *log* file). */

  norm ne = normf (ec);
  fprintf (stderr, "%d %g %g %g %g %g %g\n", 
	   LEVEL, LAPLACE, 
	   normf(un).max*sqrt(DIAMETER), 
	   ne.avg, ne.rms, ne.max,
	   ekmax);
}

#if 0
event gfsview (i += 10) {
  static FILE * fp = popen ("gfsview2D spurious.gfv", "w");
  output_gfs (fp);
}
#endif

/**
We use an adaptive mesh with a constant (maximum) resolution along the
interface. */

#if TREE
event adapt (i <= 10; i++) {
  adapt_wavelet ({c}, (double[]){0}, maxlevel = LEVEL, minlevel = 0);
}
#endif

/**
## Results

The maximum velocity converges toward machine zero for a wide range of
Laplace numbers on a timescale comparable to the viscous dissipation
timescale, as expected.

~~~gnuplot Evolution of the amplitude of the capillary currents $\max(|\mathbf{u}|)(D/\sigma)^{1/2}$ as a function of non-dimensional time $\tau=t\mu/D^2$ for the range of Laplace numbers indicated in the legend.
set xlabel 't{/Symbol m}/D^2'
set ylabel 'U(D/{/Symbol s})^{1/2}'
set logscale y
plot 'La-120-5' w l t "La=120", 'La-1200-5' w l t "La=1200", \
  'La-12000-5' w l t "La=12000"
~~~

The equilibrium shape and curvature converge toward the exact shape
and curvature at close to second-order rate.

~~~gnuplot Convergence of the error on the equilibrium shape of the droplet with resolution. The diameter is given in number of grid points.
set xlabel 'D'
set ylabel 'Shape error'
set logscale x
set xtics 2
set pointsize 1
plot [5:120]'< sort -n -k1,2 log' u (0.8*2**$1):5 w lp t "RMS", \
            '< sort -n -k1,2 log' u (0.8*2**$1):6 w lp t "Max", \
             0.2/(x*x) t "Second order"
~~~

~~~gnuplot Convergence of the relative error on the equilibrium curvature value with resolution. The diameter is given in number of grid points.
set ylabel 'Relative curvature error'
plot [5:120]'< sort -n -k1,2 log' u (0.8*2**$1):($7/2.5) w lp t "Max", \
             0.6/(x*x) t "Second order"
~~~

## See also

* [Same test with Gerris](http://gerris.dalembert.upmc.fr/gerris/tests/tests/spurious.html)
*/