# Title: Circular droplet in equilibrium
#
# Description:
#
# A circular droplet of diameter $D=0.8$ is initialised centered on
# the top-left corner of the unit box. Surface tension is imposed on
# the interface. The exact solution is given by Laplace's law: uniform
# zero velocity and a pressure jump accross the interface exactly
# balancing the surface tension force.
#
# The initial condition -- while close to the exact solution -- does
# not guarantee the exact balance between the numerical
# discretisations of surface tension and pressure gradient. However,
# these small initial perturbations generate small capillary waves
# which are progressively (on a timescale of order $D^2/\mu$) damped
# by viscosity so that the exact (to round-off error) balance is
# eventually obtained.
#
# The convergence is obtained for a wide range of Laplace numbers
# $La=\sigma\rho D/\mu^2$, as illustrated on Figure \ref{laplace}.
# Correspondingly, convergence of the curvature to a constant value is
# also obtained at all Laplace numbers as illustrated on Figure
# \ref{curvature}.
#
# Figure \ref{convergence} illustrates the convergence of the error on
# the droplet shape as a function of resolution for a Laplace number
# of 12000. Both the shape error and the relative error on the
# equilibrium curvature value illustrated on Figure
# \ref{kconvergence} show close to second-order convergence.
#
# \begin{figure}[htbp]
# \caption{\label{laplace}Evolution of the amplitude of the capillary currents
# $\max(|{\bf 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.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{laplace.eps}
# \end{center}
# \end{figure}
#
# \begin{figure}[htbp]
# \caption{\label{curvature}Evolution of the standard deviation of the
# value of the curvature along the interface as a function of
# non-dimensional time $\tau=t\mu/D^2$ for the range of Laplace
# numbers indicated in the legend.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{curvature.eps}
# \end{center}
# \end{figure}
#
# \begin{figure}[htbp]
# \caption{\label{convergence}Convergence of the error on the equilibrium shape of the
# droplet with resolution. The diameter is given in number of grid
# points.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{convergence.eps}
# \end{center}
# \end{figure}
#
# \begin{figure}[htbp]
# \caption{\label{kconvergence}Convergence of the relative error on the
# equilibrium curvature value with resolution. The diameter is given
# in number of grid points.}
# \begin{center}
# \includegraphics[width=0.8\hsize]{kconvergence.eps}
# \end{center}
# \end{figure}
#
# Author: St\'ephane Popinet
# Command: sh spurious.sh spurious.gfs 1e-11
# Version: 1.1.2
# Required files: spurious.sh convergence.ref kconvergence.ref
# Generated files: laplace.eps curvature.eps convergence.eps kconvergence.eps
#

Define CIRCLE (ellipse (-0.5,0.5,0.4,0.4))
Define MU sqrt(0.8/LAPLACE)

1 0 GfsSimulation GfsBox GfsGEdge {} {
  Time { end = TMAX }
  Refine LEVEL
  RefineSurface {return 10;} CIRCLE

  ApproxProjectionParams { tolerance = 1e-6 }
  ProjectionParams { tolerance = 1e-6 }
  AdvectionParams { scheme = none }

  VariableTracerVOF T
  VariableCurvature K T
  SourceTension T 1 K
  SourceDiffusion U MU
  SourceDiffusion V MU

  InitFraction T CIRCLE
  Init {} { Tref = T }

  AdaptGradient { istep = 1 } { cmax = 1e-6 maxlevel = LEVEL } T
  EventStop { istep = 10 } T DT

  OutputSimulation { start = end } stdout { depth = LEVEL }
  OutputScalarNorm { istep = 1 } {
    awk '{ print MU*$3/(0.8*0.8), $9*sqrt(0.8) }' > La-LAPLACE-LEVEL
  } { v = Velocity }
  OutputScalarNorm { istep = 1 } {
    awk '{ print MU*$3/(0.8*0.8), $5, $7, $9 }' > E-LAPLACE-LEVEL
  } { v = (Tref - T) }
  OutputScalarStats { istep = 1 } {
    awk '{ print MU*$3/(0.8*0.8), $5, $7, $9, $11 }' > K-LAPLACE-LEVEL
  } { v = (K - 2.50771) }
  OutputScalarNorm { istep = 1 } {
    awk '{ print MU*$3/(0.8*0.8), $5, $7, $9 }' > EK-LAPLACE-LEVEL
  } { v = (T > 0 && T < 1 ? K - 2.5 : 0) }
}
GfsBox {
  top = Boundary
  left = Boundary
}