Help with customizing boundary condition

[Michael Zhong]

Hello everyone,

I am trying to use basilisk to simulate homogeneous shear turbulence. 

Let's say shear S is on y direction and the mean velocity field is only along x direction satisfying: U_x=S*y, and the domain size is (L_x,L_y,L_z). The algorithm only simulates the fluctuation velocity field. In order to superimpose the shear effects, two key parts are necessary in the algorithm: assigning shear periodic boundary condition, and shear remapping. (The boundary conditions for x,z direction are still periodic.)

1. Assign shear periodic boundary condition 

All quantities of interest (u_x,u_y,u_z,p) should satisfy: f(x, L_y, z)=f(x-S*dt*L_y, 0, z) on top boundary; f(x, L_y, z)=f(x+S*dt*L_y, 0, z) on bottom boundary.

2. Shear remapping

After solving the fluctuation velocity field based on the shear periodic boundary condition, all those quantities should be convected based on the mean velocity field. It is done by shear remapping: u^{new}_x(x, y, z) = u_x(x-S*dt*y, y, z), u^{new}_y(x, y, z) = u_y(x-S*dt*y, y, z), u^{new}_z(x, y, z) = u_z(x-S*dt*y, y, z). 

In order to have the capability to interpolate value at x = x-S*dt*L_y or x = x-S*dt*y, I first gather the quantities across all cores through MPI and then do interpolation based on the full velocity or pressure field. 

The issue I am facing right now is that I don't know how to enforce the shear periodic boundary on top and bottom boundaries for pressure in the multigrid poisson solver. If I don't assign the top and botttom boudary conditions, it will make the output after solving the poisson equation wrong. 

In order to debug, I write two solvers in MAC scheme: one is for HIT and the other is for HST. The HIT one use the basilisk internal function periodic(bottom); but I use my customizing boundary condition for the HST one. I use MPI with 8 cores and standard basilisk Makefile for compilation. I output the face-center velocity field step by step, and I find that its maximum error in HST solver is only 1e-15 before solving poisson equation. But the maximum error increases up to 1e-3 after solving the poisson equation. Both codes are attached to this message.

I would be appreciate any suggestions. Thanks!

Appendix:

1. The full algorithm is shown as follows (ref: 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.)