Help with doing interpolation

[Michael Zhong]

Hello everyone, 

I am trying to use basilisk to simulate homogeneous shear turbulence following the paper: 

Kasbaoui, M. Houssem, et al. "An algorithm for solving the Navier–Stokes equations with shear-periodic boundary conditions and its application to homogeneously sheared turbulence." Journal of Fluid Mechanics 833 (2017): 687-716. 

I attach my code as well. There are two tricky parts: 

1. I need to assign the so-called shear periodic boundary condition by myself (line: 48-113). To make sure it satisfies such boundary condition, I need to use global interpolation to calculate the value at certain location, which might be far away from the original location. But I realize that using interpolation in foreach loop would fail if I use mpi for parallization. It seems that each core only has parts of the velocity field but I need the full velocity field to do interpolation. So my question is what would be the best way to achieve this? 

foreach_boundary (top) {

    // find the corresponding x,y,z at the bottom boundary

    // if x - shift < X0, we need to apply left,right periodic bc

    double x_interp = X0 + fmod((x - shift - X0 + L0), L0);

    double y_interp = Y0 + Delta/2;

    // u.x[0,1] = interpolate(u.x, x_interp, y_interp, z);

    u.y[0,1] = interpolate(u_temp.y, x_interp, y_interp, z);

    // u.z[0,1] = interpolate(u.z, x_interp, y_interp, z);

    uf.y[] = 0.5*(u.y[] + u.y[0,1]);

    uf.y[0,1] = interpolate(u_temp.y, x_interp, Y0 + Delta, z);

    // check interpolate value

    if (fabs(uf.y[]) > 1e6 || fabs(uf.y[0,1]) > 1e6) {

      fprintf(stderr, "Top: x,y,z,uf.y[],uf.y[0,1]:%g,%g,%g,%g,%g\n", x,y,z,uf.y[],uf.y[0,1]);

    }

  }

  foreach_boundary (bottom) {

    // find the corresponding x,y,z at the top boundary

    // if x + shift > X0 + Lx, we need to apply left,right periodic bc

    double x_edge = X0 + fmod((x + shift - X0 + L0), L0);

    double y_edge = Y0 + L0 - Delta/2; // y center at the top cell

    // u.x[0,-1] = interpolate(u.x, x_edge, y_edge, z);

    u.y[0,-1] = interpolate(u_temp.y, x_edge, y_edge, z);

    // u.z[0,-1] = interpolate(u.z, x_edge, y_edge, z);

    uf.y[] = 0.5*(u.y[] + u.y[0,-1]);

    uf.y[0,-1] = interpolate(u_temp.y, x_edge, Y0 + L0 - Delta, z);

    // check interpolate value

    if (fabs(uf.y[]) > 1e6 || fabs(uf.y[0,-1]) > 1e6) {

      fprintf(stderr, "Bottom: dt,x_edge,x,y,z,uf.y[],uf.y[0,-1],u.y[],u.y[0,-1]:%g,%g,%g,%g,%g,%g,%g,%g,%g\n", dt,x_edge,x,y,z,uf.y[],uf.y[0,-1],u.y[],u.y[0,-1]);

    }

  }

}

 My current idea is to first create several double arrays to save the interpolated value. Then use foreach_boundary to copy those value back. But I am not sure is there any better way to do this.

2. I need to do so-called shear remapping in the end of each time step (line 367-398). It is very similar to assign the shear periodic boundary condition. But this time I need to implement that for all the interior points. I also want to learn how to do the interpolation efficiently in this part.

// Apply shear remapping

//  f(x,y,z) = f(x-S*dt*y,y,z)

//  cannot apply interpolate to face vector

event shear_remap (i++, last) {

    get_cell_centered_velocity(u, uf);

   

    vector u_temp[];

    boundary ({u,u_temp,uf});

    foreach() {

      foreach_dimension() {

        u_temp.x[] = u.x[]; // Copy cell-centered velocity to temporary vector

      }

    }

    // Apply shear remapping to the cell-centered velocity

    foreach() {

        double shift = fmod(dt * SHEAR * y, L0);

        double x_interp = X0 + fmod((x - shift - X0 + L0), L0);

        u.x[] = interpolate(u_temp.x, x_interp, y, z);

        u.y[] = interpolate(u_temp.y, x_interp, y, z);

        u.z[] = interpolate(u_temp.z, x_interp, y, z);

    }

    foreach_face() {

        uf.x[] = face_value(u.x, 0);

    }

   

    apply_shear_bc(uf);

    // project to enforce incompressibility

    project(uf, p, alpha, dt, 4);

}

My current idea is to first create a double array to save all the interpolated value for each grid points and then copy them back using foreach. But I am not sure is there any better way to do this.

I would be appreciate any suggestions. Thanks!

[Michael Zhong]

Here is the code. I use:

make HST_mac.tst -f Makefile 

for compilation. The Makefile is a standard basilisk Makefile.