sandbox/Antoonvh/trefoil4.c

    Andres Castillo’s Trefoil example…

    …Using a 4th order solver and bwatch. See Andres’ page for my inspiration.

    Volumetric rendering of the negative \lambda_2 field (via surfdrive)

    #include "grid/octree.h"
    #include "nsf4t.h"
    scalar * tracers = NULL;
    #include "filaments.h"
    #include "lambda2.h"
    #include "bwatch.h"
    
    int n_seg = 1000;
    double ue = 1e-3, vis = 5e-4;
    
    int main(int argc, char ** argv) {
      if (argc > 1)
        ue = atof (argv[1]);
      if (argc > 2)
        vis = atof (argv[2]);
      foreach_dimension()
        periodic(left);
      L0 = 16;
      X0 = Y0 = Z0 = -L0/2;
      N = 1 << 5;
      a_LO = 0.2;
      const scalar visc[] = vis;
      nu = visc;
      run();
    }

    The Knot is defined as a parametric curve

    coord knot (double theta) {
      coord C;
      C.x = sin(theta) + 2.*sin(2.*theta);
      C.y = cos(theta) - 2.*cos(2.*theta);
      C.z = -sin(3.*theta) - 4;
      return C;
    }
    
    event init (t = 0) {
      refine (sq(x) + sq(y) + sq(z) < sq(4) && level < 6);
      vector omega[], psi[];
      int n_seg = 1000;
      foreach_dimension()
        psi.x.prolongation = refine_4th;
      foreach() {
        foreach_dimension()
          psi.x[] = 0;
      }
      boundary ((scalar*){psi});
      double oe = 0.01/sqrt(a_LO);
      do {
        get_vor_vector (omega, knot, 0 , 2*pi, n_seg, Lamb_Oseen);
        printf ("#cells: %ld \n", grid->tn);
      } while (adapt_wavelet ((scalar*){omega}, (double[]){oe, oe, oe}, 9).nf > (grid->tn/100));
      foreach_dimension() {
        stats so = statsf (omega.x);
        foreach()
          omega.x[] -= so.sum/so.volume;
      }
      foreach_dimension()
        poisson (psi.x, omega.x);
      vector uc[];
      foreach_dimension()
        uc.x.prolongation = refine_4th;
      foreach(){
        foreach_dimension()
          uc.x[] = -((8*(psi.z[0,1] - psi.z[0,-1]) + psi.z[0,-2] - psi.z[0,2]) -
    		 (8*(psi.y[0,0,1] - psi.y[0,0,-1]) + psi.y[0,0,-2] - psi.y[0,0,2]))/(12*Delta);
      }
      boundary ((scalar *){uc});
      vector_to_face (uc);
      project (u, p);
    }
    
    #define FUNC(x) (exp(-x) + x - 1)
    
    event mov (t += 0.1) {
      static FILE * fp = popen ("ppm2mp4 tref.mp4", "w");
      vector uc[];
      face_to_vector (uc);
      scalar l2[];
      lambda2 (uc, l2);
      foreach() 
        l2[] = l2[] < 0 ? FUNC(-l2[]): 0;
      boundary ({l2});
      watch (fov = 12, O = {18*sin(sin(t/50)), 5*sin(t/50), 20*cos(t/50)},
    	 poi = {0, 0, 2},
    	 nx = 1024, ny = (40*24));
      volume (l2, sc = .4, min = -15, max = 15,
    	  cols = true, shading = 1, mval = 1e-6);
      store (fp);
      plain();
    }
    
    event adapt (i++)
      adapt_flow (ue, 99, 1);
    
    event logger (i++) {
      fprintf (stderr, "%d %g %d %d %d %d %ld %d\n", i, t, mgp.i, 
               mgp.nrelax, mgp2.i, mgp2.nrelax, grid->tn, grid->maxdepth);
    }
    
    event stop (t = 35);