sandbox/Antoonvh/trefoil4.c
Andres Castillo’s Trefoil example…
…Using a 4th order solver and bwatch. See Andres’ page for my inspiration.
The knot unties itself
#define RKORDER 3
#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)
= atof (argv[1]);
ue if (argc > 2)
= atof (argv[2]);
vis foreach_dimension()
periodic(left);
= 16;
L0 = Y0 = Z0 = -L0/2;
X0 = 1 << 5;
N = 0.2;
a_LO const scalar visc[] = vis;
= visc;
nu run();
}
The Knot is defined as a parametric curve
knot (double theta) {
coord ;
coord C.x = sin(theta) + 2.*sin(2.*theta);
C.y = cos(theta) - 2.*cos(2.*theta);
C.z = -sin(3.*theta) - 4;
Creturn 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()
.x.prolongation = refine_4th;
psiforeach() {
foreach_dimension()
.x[] = 0;
psi}
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()
.x[] -= so.sum/so.volume;
omega}
foreach_dimension()
(psi.x, omega.x);
poisson vector uc[];
foreach_dimension()
.x.prolongation = refine_4th;
ucforeach(){
foreach_dimension()
.x[] = -((8*(psi.z[0,1] - psi.z[0,-1]) + psi.z[0,-2] - psi.z[0,2]) -
uc(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);
(u, p);
project }
#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[] < 0 ? FUNC(-l2[]): 0;
l2boundary ({l2});
watch (fov = 12, O = {18*sin(sin(t/50)), 5*sin(t/50), 20*cos(t/50)},
= {0, 0, 2},
poi = 1024, ny = (40*24));
nx volume (l2, sc = .4, min = -15, max = 15,
= true, shading = 1, mval = 1e-6);
cols 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,
.nrelax, mgp2.i, mgp2.nrelax, grid->tn, grid->maxdepth);
mgp}
event stop (t = 35);