Open boundary conditions.

[cpviolator]

Some computations (Lüscher's open temporal for example) require open boundary conditions. QUDA needs an open BC option to support this. A first pass would be to identify the open BC and insert a zero valued link as needed. An optimised version would do this, and not communicate at all in the open dimension.

[srijitpaul]

Hello @cpviolator, has there been any development in this issue?

[cpviolator]

@srijitpaul Not as yet. I have not found a use case where I have needed to implement it. Does your code require the use of OpenBC? Happy to work with you on a QUDA implementation.

[srijitpaul]

Hi Dean, yes I would require for the solver to be used on lattices with OpenBC, the list of relevant files I need to consider would be very helpful for a starting point.

[cpviolator]

@srijitpaul When I looked at this problem (and in the near future I may have to implement an open temporal boundary ) my thoughts were to have three types of boundary condition in QUDA, each of which multiply by the given number when crossing the T=0 boundary:

1. Periodic -> 1.0
2. Anti-periodic -> -1.0
3. Open -> 0.0

This would be an inefficient hack as we would still be communicating data across the T boundary despite not having to use it. Later, one could apply conditionals on the communications dependent on the BC such that if the problem is MPI split in only the T dimension, comms would be minimised.

From here I must sheepishly defer to @maddyscientist. Can you point us to a specific point in the code where the Wilson operator has BC applied?

[mathiaswagner]

But you would still have to communicate in T direction when you do not cross the T=0 boundary, so the gain would be only for nodes communicating over the T=0 <-> T=N-1 boundary. Or do I miss something here?

[cpviolator]

@mathiaswagner For sure. I was conflating this with a project that has zero temporal links!

[mathiaswagner]

We already have

  typedef enum QudaTboundary_s {
    QUDA_ANTI_PERIODIC_T = -1,
    QUDA_PERIODIC_T = 1,
    QUDA_INVALID_T_BOUNDARY = QUDA_INVALID_ENUM
  } QudaTboundary;

which could be easily extended.

A lot of the boundary stuff is handled in gauge_field_order.h. Probably a lot of the code is generic enough, but it has not been tested for open BC. It would probably also be good to add it to the reference codes for the dslashes.

[maddyscientist]

I think the easiest and fastest approach would be to just zero out the gauge field at the boundary, and run the code as is and things should just work.

@srijitpaul in general, are the open boundary conditions applied only for fermions and does this apply to the gauge action as well?

[srijitpaul]

@maddyscientist @cpviolator @mathiaswagner thanks for your replies. We have open boundary conditions for fermions and gauge action.

[maddyscientist]

Ok, that make things easier. In that case, adding open boundary conditions should just be a case of 1.) Adding an appropriate enum type to the QudaTboundary_s struct 2.) Zero out the gauge field elements on the last time slice prior to any halo exchange in GaugeField::exchangeGhost and GaugeField::exchangeExtendedGhost.

I think this alone should be sufficient. Perhaps the zeroing of the gauge field matrices on the last time slice should be done as part of GaugeField::exchangeGhost and GaugeField::exchangeExtendedGhost?

Does this sound about right? Of course, we could optimize the Dirac operators to not communicate across the temporal boundary, though I don't expect the gain to be significant. That would definitely be a secondary step.

[srijitpaul]

Ok, I have tried adding these lines in GaugeFiled::exchangeGhost with an if condition of open boundary condition, after extractGaugeGhost,

      if (this->Ncolor() == 3 && link_dir == 1) {
        if ((comm_rank()+1)%comm_dim(nDim-1) == 0) {
          qudaMemset(send_d[nDim-1], 0, ghost_face_bytes_aligned[nDim-1]);
        }
      }

The first Ncolor() condition makes sure that this condition is not used within different levels in the multigrid. The link_dir condition ensure only forward time direction gauge field are considered. The comm_rank condition ensures only the last timeslice at the boundary is considered for zeroing. For the GaugeField::exchangeGhost , I have a similar implementation,

    if (this->Ncolor() == 3) {
      if ((comm_rank()+1)%comm_dim(nDim-1) == 0) {
        qudaMemset(send_d[nDim-1], 0, 2 * ghost_face_bytes_aligned[nDim-1]);
      }
    }

Am I doing something wrong here @maddyscientist ?

[maddyscientist]

Hi @srijitpaul. Sorry for the slow reply. The code I use for checking if we are on the first / last time slice is in gauge_field_order.h:

          isFirstTimeSlice(comm_coord(3) == 0 ? true : false),
          isLastTimeSlice(comm_coord(3) == comm_dim(3) - 1 ? true : false),

Not sure if that's your issue or not, I can look more closely tomorrow.