droplet bagmode breakup

[张祥]

hello, i m a student new to basilisk,now i m simulating the droplet bagmode breakup under left flow. i have learned bagmode program from pairetti's web and  kaitaotang.and simulate successfully under low density ratio.
but i find my droplet's surface unstable when i use high density ratio,such as water 1000 and air 1.2 .i wonder how pairetti and kaitaotang made it successful to simulate in their papers. should i modify my dt? or cfl number ? and how?

above jpg is kaitaotang's result, has smooth profile under large density ratio

above is my result which has unstable surface 
and here is my settings

#define dimension 2

#define SQUARES

// #include "grid/octree.h"

#include "navier-stokes/centered.h"

#include "two-phase.h"

#include "navier-stokes/conserving.h"

#include "tension.h"

#include "view.h"

// ===========调用MPI=====================

#include <mpi.h>

#include <stdio.h>

// ============参数===================

#define RHOL 1000.

#define MUL 1.e-3

#define SIGMA 0.0728

/**Air properties*/

#define RHOG 1.2

#define MUG 1.81e-5

/**Drop diameter and stream velocity*/

#define DIAMETER 13.89e-3

#define USTREAM 8.05

int minlevel = 5;  

int MAXLEVEL = 9;  

double maxruntime = HUGE;

u.n[left]  = dirichlet(USTREAM);

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

p[left]    = neumann(0);

u.n[right] = neumann(0);

p[right]   = dirichlet(0);

// =========================================

int main (int argc, char * argv[]) {

  // 获取当前进程的编号（rank）

  int rank;

  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  // 获取总的进程数量

  int size;

  MPI_Comm_size(MPI_COMM_WORLD, &size);

  // 输出当前进程的编号和总进程数量

  printf("Hello from process %d of %d\n", rank, size);

  MPI_Barrier(MPI_COMM_WORLD);  

  L0 = 10*DIAMETER;

  X0 = -2.*DIAMETER;

  Y0 = -5.*DIAMETER;

  Z0 = -5.*DIAMETER;

  if (pid() == 0) {

    // fprintf(stderr, "Domain X: [%g, %g]\n", X0, X0 + L0);

    // fprintf(stderr, "Domain Y: [%g, %g]\n", Y0, Y0 + L0);

    // fprintf(stderr, "Domain Z: [%g, %g]\n", Z0, Z0 + L0);

    // 计算无量纲数

    double Re = RHOG * USTREAM * DIAMETER / MUG;

    double We = RHOG * USTREAM * USTREAM * DIAMETER / SIGMA;

    double Oh = MUL / sqrt(RHOL * SIGMA * DIAMETER);

    double rho_ratio = RHOL / RHOG;

    double mu_ratio  = MUL / MUG;

    // 输出

    fprintf(stderr, "\n==== Dimensionless Numbers ====\n");

    fprintf(stderr, "Reynolds number (Re):       %.2f\n", Re);

    fprintf(stderr, "Weber number (We):          %.2f\n", We);

    fprintf(stderr, "Ohnesorge number (Oh):      %.4f\n", Oh);

    fprintf(stderr, "Density ratio (rho_l/rho_g): %.2f\n", rho_ratio);

    fprintf(stderr, "Viscosity ratio (mu_l/mu_g): %.2f\n", mu_ratio);

    fprintf(stderr, "===============================\n\n");

  }

  double minDelta = L0/pow(2.0,MAXLEVEL);

  double sigmaTSLimit = sqrt(min(RHOL, RHOG)*pow(minDelta,3.0)/(SIGMA));

  double muTSLimit = min(MUL, MUG)*minDelta/SIGMA;

  DT = min(sigmaTSLimit,muTSLimit);

  // DT = 1.25e-5;

  CFL = 0.4;

  f.sigma = SIGMA;

  rho1 = RHOL;

  rho2 = RHOG;

  mu1 = MUL;

  mu2 = MUG;

  TOLERANCE = 1e-4;

  run();

 

}

event init (t = 0) {

  // if (!restore (file = "dump-0.04")) {

    refine (sq(x) + sq(y) + sq(z) - sq(0.6*DIAMETER) < 0 && sq(x) + sq(y) + sq(z) - sq(0.4*DIAMETER) > 0 && level < MAXLEVEL);

    fraction (f, sq(0.5*DIAMETER) - (sq(x) + sq(y) + sq(z)));

  // }

}

event adapt (i++) {

  double uemax = 5e-2;

  adapt_wavelet ({f,u}, (double[]){1e-3,uemax,uemax,uemax}, MAXLEVEL, minlevel);

}

event end (t=0.1) {

}

//==================输出信息=================

event logfile (i += 100 ,first) {

  if (i == 0){

    fprintf (ferr,

       "t dt mgp.i mgpf.i mgu.i perf.t perf.speed grid->tn\n");

  }

  fprintf (ferr,

     "%g %g %d %d %d"

     "%.2e %.2e %ld \n",

     t, dt, mgp.i, mgpf.i, mgu.i,

     perf.t, perf.speed, grid->tn);

  fflush (ferr);

}

//===================输出图像=========================

//====================================

int frame = 0;

event movies ( t += 1.e-3 )

{

  view (fov = 12, width = 1920, height = 1080,tx = -0.4);

  clear();

  draw_vof ("f", filled = 1, lw = 1, lc = {1, 1, 1}, fc = {0.6, 0.8, 1.0});

  cells();

  char name[256];

 

  sprintf(name, "ppm/snapshot-%d.ppm", frame++);

  save(name);

}

[张祥]

under high density ratio,my surface unstable again. i have used the dt in pairetti's programe,eg.

  double minDelta = L0/pow(2.0,MAXLEVEL);

  double sigmaTSLimit = sqrt(min(RHOL, RHOG)*pow(minDelta,3.0)/(SIGMA));

  double muTSLimit = min(MUL, MUG)*minDelta/SIGMA;

  DT = min(sigmaTSLimit,muTSLimit);

what should i do to solve this problem?