sandbox/hasansh/fast-capwave/fast-capwave-density.c
Capillary wave for fluids with different densities using new poisson solver
We use a constant-resolution grid, the Navier–Stokes solver with VOF interface tracking and surface tension. Since we’re using the new fast poisson solver we include fast-navierstokes.h instead of navier-stokes/centered.h.
#include "grid/multigrid.h"
#include "fast-navierstokes.h"
#include "two-phase.h"
#include "tension.h"
#include "prosperetti-density.h"
We make sure that the boundary conditions for the face-centered velocity field are consistent with the centered velocity field (this affects the advection term).
uf.n[left] = 0.;
uf.n[right] = 0.;
uf.n[top] = 0.;
uf.n[bottom] = 0.;
We will store the accumulated error in se and the number of samples in ne.
double se = 0; int ne = 0;
int main() {
The domain is 2x2 to minimise finite-size effects. The surface tension is one and the viscosity is constant.
L0 = 2.;
Y0 = -L0/2.;
f.sigma = 1;
TOLERANCE = 1e-6;
mu1 = 0.0182571749236;
mu2 = 0.0182571749236;
The density of fluid 1 is one and the density ratio is 10.
rho1 = 1;
rho2 = 0.1;
\rho_0 is defined as minimum of \rho_1 and \rho_2.
rho0 = min(rho1, rho2);
We change the methods to compare the results.
system ("rm -f convergence-*");
for (int iter = unsplit; iter <= fastpstar; iter++) {
We vary the resolution to check for convergence.
for (N = 16; N <= 128; N *= 2) {
se = ne = 0;
run();
}
poisson_solver = iter + 1;
}
}
The initial condition is a small amplitude plane wave of wavelength unity.
By default tracers are defined at t-\Delta t/2. We use the first keyword to move VOF advection before the amplitude output i.e. at t+\Delta/2. This improves the results.
event vof (i++, first);
We output the amplitude at times matching exactly those in the reference file.
To get an accurate amplitude, we reconstruct the height function field and take the corresponding maximum.
vector h[];
heights (f, h);
double max = - HUGE;;
foreach()
if (f[] > 0 && f[] < 1) {
double yi = y + height(h.y[])*Delta;
if (yi > max)
max = yi;
}
We output the corresponding evolution in a file indexed with the method poisson_solver and the number of grid points N.
char name[80];
sprintf (name, "wave-%i-%d", poisson_solver ,N);
static FILE * fp = fopen (name, "w");
fprintf (fp, "%g %g\n", t*15.016663878457, max);
fflush (fp);
To compute the RMS error, we get data from the reference file prosperetti-density.h and add the difference to the accumulated error.
se += sq(max - prosperetti[ne][1]); ne++;
}
At the end of the simulation, we output on a file, indexed with the method poisson_solver, the resolution (number of grid points per wavelength) and the relative RMS error.
event error (t = end) {
char name[80];
sprintf (name, "convergence-%i", poisson_solver);
FILE * fp = fopen (name, "a");
if (fp != NULL) {
fprintf (fp, "%g %g\n", N/L0, sqrt(se/ne)/0.01);
fflush (fp);
}
else {
FILE * fp = fopen (name, "w");
fprintf (fp, "%g %g\n", N/L0, sqrt(se/ne)/0.01);
fflush (fp);
}
}
Results
set xlabel 'tau'
set ylabel 'Relative amplitude'
plot '../prosperetti-density.h' u 2:4 w l t "Prosperetti", \
'wave-0-128' every 10 w p t "Basilisk",\
'wave-1-128' every 10 w p t "FastPn",\
'wave-2-128' every 10 w p t "FastP*"
set xlabel 'Number of grid points'
set ylabel 'Relative RMS error'
set logscale y
set logscale x 2
set grid
plot [5:200][1e-4:1] 2./x**2 t "Second order",\
'convergence-0' t "Basilisk" w lp, \
'convergence-1' t "FastPn" w lp, \
'convergence-2' t "FastP*" w lp