src/examples/hele-shaw.c

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    
    #include "hele-shaw.h"
    
    #define MAXLEVEL 8
    
    double mu1 = 2.5, mu2 = 1e-3, k = 0.1;
    
    scalar f[];
    scalar * tracers = {f};
    
    p[left]  = dirichlet(1e-3);
    p[right] = dirichlet(0);
    f[left]  = 1.;
    
    int main()
    {
      size (7.5e-2);
      init_grid (1 << MAXLEVEL);
      run();
    }
    
    event init (i = 0) {
      foreach()
        f[] = (x < 1e-3)*(1. - 1e-3*(1. + noise()));
    }
    
    event logfile (i++)
    {
      stats s = statsf (f);
      fprintf (stderr, "%d %g %d %g %g %g\n", 
    	   i, t, mgp.i, s.sum, s.min, s.max);
    }
    
    event movies (t += 0.2)
    {
      output_ppm (f, min = 0, max = 1, file = "f.mp4");
    
      output_ppm (u.y, file = "v.mp4");
    
      scalar l[];
      foreach()
        l[] = level;
      output_ppm (l, min = 0, max = MAXLEVEL, linear = false, file = "level.mp4");
    }
    
    event output (t += 10; t <= 50)
    {
      static int nf = 0;
      printf ("file: field-%d\n", nf++);
      output_field ({p,f,u}, stdout, N);
    }
    
    event coefficients (i++)
    {
      foreach_face() {
        double ff = (f[] + f[-1])/2.;
        beta.x[] = - k/(mu1 + clamp(ff,0,1)*(mu2 - mu1));
      }
    }
    
    #if TREE
    event adapt (i++) {
      double tolerance = 0.02;
      adapt_wavelet ({f}, &tolerance, MAXLEVEL);
      event ("coefficients");
    }
    #endif