sandbox/Antoonvh/moist_bubble.c

    A rising moist bubble

    A moist bubble is placed in the atmosphere and rises because water vapor is lighter than air. It reaches condensation levels soon to form a cloud. The setup is inspired by Grabowski and Clark (1991) and the setup of B van Stratum.

    //#include "grid/octree.h"
    #include "navier-stokes/centered.h"
    #include "thermo.h"
    
    int maxlevel = 9;
    
    int main(int argc, char * argv[]) {
      if (argc > 1)
        maxlevel = atoi (argv[1]);
      const face vector K[] = {0.5, 0.5};
      mu = K;
      L0 = 3600;
      X0 = Z0 = -L0/2;
      periodic (left);
    #if (dimension == 3)
      periodic (front);
    #endif
      N = 64;
      run();
    }
    
    #define RADIUS (sqrt(sq(x) + sq(y - 800.) + sq (z)))
    event init (t = 0) {
      T_ref = 283.;
      P0 = 1e5;
      double dlnthetadz = 1.0e-5;
      refine (RADIUS < 350. && level < maxlevel);
      foreach() 
        thl[] = exp(log(283.) + dlnthetadz*y) + 0.05*noise();

    q_t may be initialized via the relative humidity (using the QSAT macro). However it requires an iterative method to find the corresponding hydrostatic pressure profile. The procedure is started after we have taken a first guess.

      scalar p_temp[];
      int j = 0, j_max = 20;
      set_pres (guess = true);
      do {
        foreach() {
          double r = RADIUS;
          if (r < 200)
    	qt[] = QSAT;
          else if (r < 300) 
    	qt[] = QSAT * (0.2 + 0.8 * sq(cos(pi*(r - 200.)/200.)));
          else 
    	qt[] = QSAT * 0.2;
        }
        set_pres ();
        j++;
      } while (change (pres, p_temp) > 0.01 && j < j_max);
      if (j == j_max)
        fprintf (stderr,
    	     "Static Pressure not found with %d iterations\n", j);
    }
    
    #include "diffusion.h"
    event tracer_diffusion (i++) {
      for (scalar s in tracers)
        diffusion (s, dt, mu);
    }
    
    event logger (i += 5) 
      fprintf(stdout, "%g %d\n", t, i);
    
    event adapt (i++) 
      adapt_wavelet ((scalar*){qt, u},
    		 (double[]){0.00025, 0.05, 0.05, 0.05}, maxlevel, 5);
    
    event stop (t = 1000);

    Output and Results

    We output movies. In 3D there there are a bunch of dumps for post processing.

    2D: The cloud

    and the 3D results, via the dumps to vtk and Paraview.

    event movie (t += 5) {
      scalar qc[], lev[];
      foreach() {
        lev[] = level;
        qc[] = QC;
      }
      output_ppm (lev, file = "level.mp4", n = 512, min = 5, max = 10);
      output_ppm (qc, file = "qc.mp4", n  = 512,
    	      map = cloud, min = 0, max = 0.0002);
      output_ppm (qt, file = "qt.mp4", n  = 512);
    }
    
    #if (dimension == 3)
    event dumps (t += 2.5) {
      scalar qc[];
      foreach()
        qc[] = QC;
      boundary ({qc});
      char fname[99];
      sprintf(fname, "dump%g", t);
      dump (fname);
    }
    #endif

    Reference

    Grabowski, Wojciech W., and Terry L. Clark. “Cloud–environment interface instability: Rising thermal calculations in two spatial dimensions.” Journal of the Atmospheric Sciences 48.4 (1991): 527-546.