adapt_wavelet_bug.c

Rising bug

Rising bug

This is a slight modification done on the rising bubble example. When compiling this example with the -catch flag, the simulation crashes as soon as adapt_wavelet is called. gdb reveals that the crash happens at:
> Program received signal SIGFPE, Arithmetic exception.
> 0x000000000047a2f6 in _boundary1 (point=…, neighbor=…, _s=…)
> at /home/gerris/basilisk/src/navier-stokes/centered.h:79
> 79 p[left] = neumann(-a.x[]*fm.x[]/alpha.x[]);
An inspection of the fields shows that a.x is not defined near the bubble (ie in the refined region).

#if 0
# include "axi.h"
#endif
#include "navier-stokes/centered.h"
#include "vof.h"
#include "tension.h"

Hysing et al. consider two cases (1 and 2), with the densities, dynamic viscosities and surface tension of fluid 1 and 2 given below.

#define rho1 1000.
#define mu1 10.

#if 1
# define rho2 100.
# define mu2 1.
# define SIGMA 24.5
#else
# define rho2 1.
# define mu2 0.1
# define SIGMA 1.96
#endif

#define LEVEL 8
#define MAXLEVEL 10

The interface between the two-phases will be tracked with the volume fraction field f. We allocate the fields for the variable density and viscosity (alphav, alphacv and muv respectively).

scalar f[];
scalar * interfaces = {f};
face vector alphav[];
scalar rhov[];
face vector muv[];

The boundary conditions are slip lateral walls and no-slip on the right and left walls.

u.t[right]  = dirichlet(0);
uf.t[right] = dirichlet(0);
u.t[left]   = dirichlet(0);
uf.t[left]  = dirichlet(0);

uf.n[top]    = 0;
uf.n[bottom] = 0;
p[top]       = neumann(0);
p[bottom]    = neumann(0);

int main() {

The domain will span [0:2]\times[0:0.5] and will be resolved with 256\times 64 grid points.

  size (2);
  init_grid (1 << LEVEL);

The density and viscosity are defined by the variable fields we allocated above. We also set the surface tension for interface f. We reduce the tolerance on the Poisson and viscous solvers to improve the accuracy.

  alpha = alphav;
  rho = rhov;
  mu = muv;
  f.sigma = SIGMA;
  TOLERANCE = 1e-4;

  run();
}

event init (t = 0) {

The bubble is centered on (0.5,0) and has a radius of 0.25.

  vertex scalar phi[];
  foreach_vertex()
    phi[] = sq(x - 0.5) + sq(y) - sq(0.25);
  fractions (phi, f);
}

We add the acceleration of gravity.

event acceleration (i++) {
  face vector av = a;
  foreach_face(x)
    av.x[] -= 0.98;
}

The density and viscosity are defined using the arithmetic average.

#define rho(f) ((f)*rho1 + (1. - (f))*rho2)
#define mu(f)  ((f)*mu1 + (1. - (f))*mu2)

event properties (i++) {
  foreach_face() {
    double ff = (f[] + f[-1])/2.;
    alphav.x[] = fm.x[]/rho(ff);
    muv.x[] = fm.x[]*mu(ff);
  }
  foreach()
    rhov[] = cm[]*rho(f[]);
}

event adapt (i = 10; i++) {
  adapt_wavelet ({f,u.x,u.y}, (double[]){1e-2,5e-3,5e-3}, MAXLEVEL);
  event ("properties");
}

At t=3 we output the shape of the bubble.

event interface (t = 0.2) {
  output_facets (f, stderr);
}