Extended gauge field routines do not handle boundary conditions correctly

[maddyscientist]

The extended gauge copy and ghost extraction routines do not presently handle the time boundary conditions correctly. When extended, the boundary condition lies on the the R[3]-1 on the first time slice of nodes and on X[3]-R[3]-1 on the last time slice of nodes (where X[] is the local extended dimensions and R[] is the extended radii dimensions).

I think to fix this, we need a modified variant of the timeBoundary functor used by the gauge field order structs. This requires a couple of additions to the GaugeField class such that any gauge field that is extended knows it is extended and it knows what the depth of the extended region is (the R[] array).

  /**
     timeBoundary variant for extended gauge field
     @param idx extended field linear index
     @param X the extended gauge field dimensions
     @param R the radii dimenions of the extended region
     @param tBoundary the boundary condition
     @param isFirstTimeSlice if we're on the first time slice of nodes
     @param isLastTimeSlide if we're on the last time slice of nodes
  */
  template <typename Float>
    __device__ __host__ inline Float timeBoundary(int idx, const int X[QUDA_MAX_DIM], const int R[QUDA_MAX_DIM], 
                          QudaTboundary tBoundary, bool isFirstTimeSlice, bool isLastTimeSlice) {
    if ( idx >= (R[3]-1)*X[0]*X[1]*X[2]/2 && idx < R[3]*X[0]*X[1]*X[2]/2 ) {
      // the boundary condition is on the R[3]-1 time slice
      return isFirstTimeSlice ? tBoundary : 1.0;
    } else if ( idx >= (X[3]-R[3]-1)*X[0]*X[1]*X[2]/2 && idx < (X[3]-R[3])*X[0]*X[1]*X[2]/2 ) { 
      // the boundary condition lies on the X[3]-R[3]-1 time slice
      return isLastTimeSlice ? tBoundary : 1.0;
    } else {
      return 1.0;
    }
  }

[AlexVaq]

Actually adding the R[] array to the

Alex

El 5/11/2014, a las 9:17, mikeaclark notifications@github.com escribió:

The extended gauge copy and ghost extraction routines do not presently handle the time boundary conditions correctly. When extended, the boundary condition lies on the the R[3]-1 on the first time slice of nodes and on X[3]-R[3]-1 on the last time slice of nodes (where X[] is the local extended dimensions and R[] is the extended radii dimensions).

I think to fix this, we need a modified variant of the timeBoundary functor used by the gauge field order structs. This requires a couple of additions to the GaugeField class such that any gauge field that is extended knows it is extended and it knows what the depth of the extended region is (the R[] array).

/ timeBoundary variant for extended gauge field @param idx extended field linear index @param X the extended gauge field dimensions @param R the radii dimenions of the extended region @param tBoundary the boundary condition @param isFirstTimeSlice if we're on the first time slice of nodes @param isLastTimeSlide if we're on the last time slice of nodes _/ template device host* inline Float timeBoundary(int idx, const int X[QUDA_MAX_DIM], const int R[QUDA_MAX_DIM], QudaTboundary tBoundary, bool isFirstTimeSlice, bool isLastTimeSlice) { if ( idx >= (R[3]-1)_X[0]_X[1]_X[2]/2 && idx < R[3]_X[0]_X[1]_X[2]/2 ) { // the boundary condition is on the R[3]-1 time slice return isFirstTimeSlice ? tBoundary : 1.0; } else if ( idx >= (X[3]-R[3]-1)_X[0]_X[1]_X[2]/2 && idx < (X[3]-R[3])_X[0]_X[1]X[2]/2 ) { // the boundary condition lies on the X[3]-R[3]-1 time slice return isLastTimeSlice ? tBoundary : 1.0; } else { return 1.0; } }

— Reply to this email directly or view it on GitHub https://github.com/lattice/quda/issues/173.

[AlexVaq]

Actually adding the R[] array to the GaugeField class might be enough. One can just set it to zero when the fields are not extended, and use the extended formulae for computing the boundaries, the sizes, etc. everywhere. The final result should remain correct for extended and non-extended gauge fields.

Alex

El 5/11/2014, a las 9:17, mikeaclark notifications@github.com escribió:

The extended gauge copy and ghost extraction routines do not presently handle the time boundary conditions correctly. When extended, the boundary condition lies on the the R[3]-1 on the first time slice of nodes and on X[3]-R[3]-1 on the last time slice of nodes (where X[] is the local extended dimensions and R[] is the extended radii dimensions).

I think to fix this, we need a modified variant of the timeBoundary functor used by the gauge field order structs. This requires a couple of additions to the GaugeField class such that any gauge field that is extended knows it is extended and it knows what the depth of the extended region is (the R[] array).

/ timeBoundary variant for extended gauge field @param idx extended field linear index @param X the extended gauge field dimensions @param R the radii dimenions of the extended region @param tBoundary the boundary condition @param isFirstTimeSlice if we're on the first time slice of nodes @param isLastTimeSlide if we're on the last time slice of nodes _/ template device host* inline Float timeBoundary(int idx, const int X[QUDA_MAX_DIM], const int R[QUDA_MAX_DIM], QudaTboundary tBoundary, bool isFirstTimeSlice, bool isLastTimeSlice) { if ( idx >= (R[3]-1)_X[0]_X[1]_X[2]/2 && idx < R[3]_X[0]_X[1]_X[2]/2 ) { // the boundary condition is on the R[3]-1 time slice return isFirstTimeSlice ? tBoundary : 1.0; } else if ( idx >= (X[3]-R[3]-1)_X[0]_X[1]_X[2]/2 && idx < (X[3]-R[3])_X[0]_X[1]X[2]/2 ) { // the boundary condition lies on the X[3]-R[3]-1 time slice return isLastTimeSlice ? tBoundary : 1.0; } else { return 1.0; } }

— Reply to this email directly or view it on GitHub https://github.com/lattice/quda/issues/173.

[maddyscientist]

I believe this issue is now fixed (as of commit a26d948b706eaa3583685cf5882789a9f0d4740a).