Question: connectivity with shear periodic boundary conditions

[stichou]

I am using p4est with the objective of building upon it a finite volume solver.
For one particular case I need to work in a shearing box of dimensions $L_x$ and $L_y$.
Specifically, this box is periodic in $y$ direction and shear-periodic in $x$ direction:
$f(x + L_x, y+ \delta y(t) ) = f(x, y) $
$f(x, y + L_y) = f(x,y)$
where $\delta y (t)= s L_x \mathrm{mod}(t, \frac{L_y}{s L_x})$ where $s$ is the shear rate assumed constant.

In other words, the $x$ boundary is shifted at each time step in the $y$ direction and eventually comes back to a fully periodic configuration when $t= n \frac{L_y}{s L_x}$ ($n \in \mathbf{N}$).

Problem

For satisfying above BC, I thought to modify the connectivity at each time step but this task seems highly inefficient and also according to https://www.p4est.org/tutorial-connectivity.html the connectivity cannot be modified during the simulation.

Possible solution

I thought that a simple solution would be to use a brick connectivity only periodic in the $y$ direction. Then I thought of extending the domain in the $x$ direction with a ghost layer and update the values according to the sheared transformation. Do you think that this is a good idea or do you have better recommendations ?

Thank you for your time.

[pkestene]

Hello @stichou,

As you state in your problem section above, I don't think you can achieve what you want to do by modifying dynamically the p4est connectivity. You need to think of the connectivity as a regular (conform) coarse mesh.

Nevertheless, you could design a p4est connectivity, a variant of the p4est brick connectivity to have shifted periodicity in one direction, but the shift would a integer multiple of the number of trees in that direction. I guess the shift rate is data dependent (and also may change in time).

I the past I implemented shearing box BC in a code using a regular mesh (no AMR) using a ghost layer for transfering data from one side to the shifted other side. The shift was time dependent, and data need to be interpolated (or somehow adapted) on the shifted other side (because shift was not a integer multiple of the cell mesh size).

Here, one additional technical point to consider, assuming you use a brick connectivity (which is only periodic in y-direction), when you perform AMR cycle (refine, coarsen, 2:1 balance),

• you will probably need to first design an algorithm which determines for each MPI process touching the left x-axis border, compute the number of quadrants at the border, and compute with which other MPI processes touching the right x-axis border it should sends data to.
• you probably need this algorithm to communicate user data between the two sides, but also to communicate refine criterium data, to make sure the local mesh sizes roughly match on both sides of your shearing box (roughly match means here, to ensure there is at most one AMR level of difference between both side of shear border).

[cburstedde]

This all makes sense, and thanks to @pkestene. If the x-periodicity is not built into the connectivity, you'd be designing a numerical update scheme yourself. There are at least three ways to do this:

1. use the ghost layer to search for neighbor quadrants and their owners.
2. use p4est_search_local and p4est_search_partition to find the surrounding element for a point.
3. construct an lnodes structure for face-only connectivity, degree = -1, and identify distinct bt geometrically matching lnodes to be the same DOF, manually.

Good luck and please let us know if we can help.

[stichou]

@pkestene and @cburstedde, thank you for your clear answers and insights.

@pkestene, thank you for sharing your experience with sheared periodic boundary conditions (BC). Indeed, I was exactly wondering how to transfer data when mesh elements are not exactly aligned.

There is still something that is not clear to me, and I take the opportunity to ask for clarification here. For the sake of clarity, let's assume that we are working with more than one processor. I understand that the purpose of p4est_ghost* functions is to create locally a layer of quadrants that are adjacent to the domain and belong to neighboring processors. These quadrants are then synchronized with the data belonging to these neighboring processors.

When the domain is periodic in all directions, this is simple since every local process will have a ghost halo. However, if one direction is not periodic (for example, the $x$ direction in the brick connectivity), there should be a situation where some quadrants are touching the domain boundary. So, I imagine that p4est_ghost_new will not create a halo at the domain boundary. My questions are:

How can I extend the domain or create a ghost layer at the domain boundary? (I'd be interested if you had an example, but I haven't seen it in the ones you propose, and in p4est_step3.c, the boundary conditions are periodic in both directions.)

I anticipate that if the construction of this ghost layer is possible at the domain boundary, a certain fixed level of refinement would be imposed. Would this be possible, or will the ghost layer necessarily inherit the refinement of the quadrants belonging to the domain and located at the boundary?

Thank you in advance!

[pkestene]

@stichou

Just a few general remarks about the tree connectivity and ghost to better understand the difference between periodic and non-periodic and the implications about how to handle quadrant user data, and finally implement shear-border periodicity.

• you mention p4est_step3 example; its is a good starting point; one important point to understand, is that, for simplicity this example uses the pointer q.p->user_data to store actual user application data. Doing so simplifies the use of p4est, because, the quadrant mesh metadata and the quadrant user application data are stored in the same data structure. It means for example that a cal to p4est_ghost_exchange_data will transfer quadrant mesh metadata as well as application user data. If you haven't done so, I encourage you to play with this example. Change for example the connectivity to be unitsquare (non-periodic); the code will crash because of the assertion P4EST_ASSERT (sides->elem_count == 2); in all the p4est_iterate call using a face-oriented callback. All the quadrants touching external border have one face or two face (the corner quadrants) for which sides->elem_count is one (no neighbor across that face). So you need to implement a special treatment for those faces, this is where the "maths" of the non-periodic border conditions need to be implement. For usual border conditions like e.g. "wall" or "absorbing border", implementing the border conditions into your numerical scheme is often easy as it only depends on local data, even often only on the current quadrant data (you usually don't need a ghost layer to do that). In case of a shear-periodic border, it is more complex because you need to implement this specific ghost layer.
• p4est_ghost* routines are designed to add a layer of quadrants all around a MPI subdomain. If you chose a connectivity (let say brick connectivity for simplicity) that is periodic, p4est_ghost* will also add a layer of quadrants all around the outer domain (and take of the 2:1 balance constraint). In other words, when you declare the connectivity to be periodic, there is no difference between dealing with a quadrant touching a MPI sub-domain border, or touching the outer domain border. Everything is taken care of by the p4est_ghost* routines, every face will have 2 sides.
• when you declare a outer border to be non-periodic, there won't be any ghost neighbor quadrants for regular quadrant touching the outer non-periodic border. So you need to do something special here. This brings us back to my first answer, and also Carsten's answer. You need a way (algorithm) to identify which quadrants are touching the outer border; among those quadrants identify with which MPI process you need to send data to, to receive data from; perform the user data gathering, perform the actual MPI send/receive, and then perform your numerical scheme. You need to implement that for the shearing box border condition; in each MPI process, design a data structure to hold data for these shear-ghost quadrants.
• you should probably decouple the different technical difficulties here, design simple test cases. First consider the case of uniform mesh (no AMR, level min = level max). Across a shear border, a given quadrant will have (from the geometrical point of view) two "shear-ghost" neighbors that partly touch this quadrant, then you'll probably need to perform some kind of user data interpolation.
• when you activate AMR, you surely want that, across a shear border, a given quadrant sees at most one AMR level difference (as for any other quadrant, thanks to the 2:1 balance algorithm). To do that you will probably need to MPI transfer additional data, so that the refine criterion computed on both side of the shear border are matching across a shear border; this is probably sufficient to ensure the 2:1 balance.

All in all, surely it will take some efforts but it's an interesting problem.

[stichou]

@pkestene

Thank you again for all the advice and detailed explanations. It is much clearer now.

I would like to add that my initial idea of adding an extra ghost layer for outer non-periodic boundaries came from the desire to have an identical treatment in the callback function for every face, regardless of whether it touches the physical boundary or not, when passed to p4est_iterate. However, you are correct that I only need local quadrant data, and at this stage of development, that will be the simpler solution.

Thank you again, and I hope to come back in a few weeks with a working example that includes shear periodic boundary conditions.

[stichou]

@pkestene @cburstedde

I am experimenting with the boundary conditions as you suggested and would appreciate confirmation on the convention used for the sides of faces. For this example, I am using a brick connectivity that is periodic in all directions, and I am only interested in the y direction. I assume there is no refinement. I implemented a simple callback function that iterates over faces perpendicular to the y direction and populates adjacent quads with -5 on the left side (side 0) and 5 on the right side (side 1). The function looks like this:

void flux(p4est_iter_face_info_t * info, void *user_data)
{
  p4est_iter_face_side_t *side[2];
  sc_array_t                     *sides = &(info->sides);
  data_t                            *ghost_data = (data_t *)user_data;
  data_t                            *data[2];
  context_t                       *ctx = (context_t *)p4est->user_pointer;
  int                                  i, face_dir;

  side[0] = p4est_iter_fside_array_index_int(sides, 0);
  side[1] = p4est_iter_fside_array_index_int(sides, 1);
  face_dir = side[0]->face / 2; /* 0 == x, 1 == y, 2 == z */

  if (1 != face_dir)
  {
      return;
  } 

  for (i = 0; i < 2; i++)
  {
      if (side[i]->is.full.is_ghost)
      {
          data[i] = &ghost_data[side[i]->is.full.quadid];
      } else {
          data[i] = (data_t *)side[i]->is.full.quad->p.user_data;
      }
  }

  data[0]->rho=-5;
  data[1]->rho=+5;

For interior faces, this works as expected. Here is an example for faces 0 and 7:

However, when the face touches an outer border, the situation reverses: side 0 becomes the right face, and side 1 becomes the left face. This is depicted in the figure below for face 13:

I noticed that a difference between interior faces and outer faces lies in the tree_boundary variable, which is 0 for the former and P4EST_CONNECT_FACE for the latter.

I was wondering if this reversal on the "left" and "right" quadrants adjacent to an outer face is to be expected every time that a border is periodic? Also, to be sure, can this reversal be detected just by checking that info->tree_boundary != 0?

[cburstedde]

Sides[0] contains the lowest z-order quadrant of all quadrants contained in sides[0] and sides[1]. It makes sense that the quadrant in the lowest numbered tree is in sides[0]. This is independent of the tree boundary status.

Note that in your program you're writing into every local quadrant multiple times: once for each of its faces. The order of these writes is undefined.