Final fix to the vector potential consistancy in serial, to do the

same as with MPI

Author: glesur

src/fluid/constrainedTransport/constrainedTransport.hpp

@@ -102,6 +102,9 @@ class ConstrainedTransport {
   void EvolveVectorPotential(real, IdefixArray4D<real> &);
   void ComputeMagFieldFromA(IdefixArray4D<real> &Vein, IdefixArray4D<real> &Vsout);
   void EnforceVectorPotentialBoundary(IdefixArray4D<real> &Vein); // Enforce BCs on A
+  void EnforceEMFBoundaryPeriodic(IdefixArray3D<real> ex,
+                                  IdefixArray3D<real> ey,
+                                  IdefixArray3D<real> ez);
 
 #ifdef WITH_MPI
   // Exchange surface EMFs to remove interprocess round off errors

src/fluid/constrainedTransport/enforceEMFBoundary.hpp

@@ -42,70 +42,62 @@ void ConstrainedTransport<Phys>::EnforceEMFBoundary() {
     #endif
   #endif
 
-  IdefixArray3D<real> ex = this->ex;
-  IdefixArray3D<real> ey = this->ey;
-  IdefixArray3D<real> ez = this->ez;
-
   // Enforce specific EMF regularisation
   for(int dir=0 ; dir < DIMENSIONS ; dir++ ) {
     if(data->lbound[dir] == shearingbox || data->rbound[dir] == shearingbox) {
       SymmetrizeEMFShearingBox();
     }
-    #ifdef ENFORCE_EMF_CONSISTENCY
-      if(data->lbound[dir] == periodic && data->rbound[dir] == periodic) {
-        // If domain decomposed, periodicity is already enforced by ExchangeAll
-        if(data->mygrid->nproc[dir] == 1) {
-          int ioffset = (dir == IDIR) ? data->np_int[IDIR] : 0;
-          int joffset = (dir == JDIR) ? data->np_int[JDIR] : 0;
-          int koffset = (dir == KDIR) ? data->np_int[KDIR] : 0;
-
-          int ibeg = (dir == IDIR) ? data->beg[IDIR] : 0;
-          int iend = (dir == IDIR) ? data->beg[IDIR]+1 : data->np_tot[IDIR];
-          int jbeg = (dir == JDIR) ? data->beg[JDIR] : 0;
-          int jend = (dir == JDIR) ? data->beg[JDIR]+1 : data->np_tot[JDIR];
-          int kbeg = (dir == KDIR) ? data->beg[KDIR] : 0;
-          int kend = (dir == KDIR) ? data->beg[KDIR]+1 : data->np_tot[KDIR];
-          idefix_for("BoundaryEMFPeriodic",kbeg,kend,jbeg,jend,ibeg,iend,
-            KOKKOS_LAMBDA (int k, int j, int i) {
-              real em;
-
-              if(dir==IDIR) {
-                em = HALF_F*(ez(k,j,i)+ez(k,j,i+ioffset));
-                ez(k,j,i) = em;
-                ez(k,j,i+ioffset) = em;
-
-                #if DIMENSIONS == 3
-                em = HALF_F*(ey(k,j,i)+ey(k,j,i+ioffset));
-                ey(k,j,i) = em;
-                ey(k,j,i+ioffset) = em;
-                #endif
-              }
-
-              if(dir==JDIR) {
-                em = HALF_F*(ez(k,j,i)+ez(k,j+joffset,i));
-                ez(k,j,i) = em;
-                ez(k,j+joffset,i) = em;
-
-                #if DIMENSIONS == 3
-                em = HALF_F*(ex(k,j,i)+ex(k,j+joffset,i));
-                ex(k,j,i) = em;
-                ex(k,j+joffset,i) = em;
-                #endif
-              }
-
-              if(dir==KDIR) {
-                em = HALF_F*(ex(k,j,i)+ex(k+koffset,j,i));
-                ex(k,j,i) = em;
-                ex(k+koffset,j,i) = em;
-
-                em = HALF_F*(ey(k,j,i)+ey(k+koffset,j,i));
-                ey(k,j,i) = em;
-                ey(k+koffset,j,i) = em;
-              }
-            });
-        }
+  }
+  #ifdef ENFORCE_EMF_CONSISTENCY
+  EnforceEMFBoundaryPeriodic(ex,ey,ez);
+  #endif //ENFORCE_EMF_CONSISTENCY
+#endif // MHD==YES
+  idfx::popRegion();
+}
+
+template<typename Phys>
+void ConstrainedTransport<Phys>::EnforceEMFBoundaryPeriodic(IdefixArray3D<real> ex,
+                                                            IdefixArray3D<real> ey,
+                                                            IdefixArray3D<real> ez) {
+  idfx::pushRegion("Emf::EnforceEMFBoundaryPeriodic");
+  #if MHD == YES
+  for(int dir=0 ; dir < DIMENSIONS ; dir++ ) {
+    if(data->lbound[dir] == periodic && data->rbound[dir] == periodic) {
+      // If domain decomposed, periodicity is already enforced by ExchangeAll
+      if(data->mygrid->nproc[dir] == 1) {
+        int ioffset = (dir == IDIR) ? data->np_int[IDIR] : 0;
+        int joffset = (dir == JDIR) ? data->np_int[JDIR] : 0;
+        int koffset = (dir == KDIR) ? data->np_int[KDIR] : 0;
+
+        int ibeg = (dir == IDIR) ? data->beg[IDIR] : 0;
+        int iend = (dir == IDIR) ? data->beg[IDIR]+1 : data->np_tot[IDIR];
+        int jbeg = (dir == JDIR) ? data->beg[JDIR] : 0;
+        int jend = (dir == JDIR) ? data->beg[JDIR]+1 : data->np_tot[JDIR];
+        int kbeg = (dir == KDIR) ? data->beg[KDIR] : 0;
+        int kend = (dir == KDIR) ? data->beg[KDIR]+1 : data->np_tot[KDIR];
+        idefix_for("BoundaryEMFPeriodic",kbeg,kend,jbeg,jend,ibeg,iend,
+          KOKKOS_LAMBDA (int k, int j, int i) {
+            if(dir==IDIR) {
+              ez(k,j,i+ioffset) = ez(k,j,i);
+              #if DIMENSIONS == 3
+              ey(k,j,i+ioffset) = ey(k,j,i);
+              #endif
+            }
+
+            if(dir==JDIR) {
+              ez(k,j+joffset,i) = ez(k,j,i);
+              #if DIMENSIONS == 3
+              ex(k,j+joffset,i) = ex(k,j,i);
+              #endif
+            }
+
+            if(dir==KDIR) {
+              ex(k+koffset,j,i) = ex(k,j,i);
+              ey(k+koffset,j,i) = ey(k,j,i);
+            }
+          });
       }
-    #endif //ENFORCE_EMF_CONSISTENCY
+    }
   }
 #endif // MHD==YES
   idfx::popRegion();

src/fluid/constrainedTransport/enforceVectorPotentialBoundary.hpp

@@ -28,6 +28,7 @@ void ConstrainedTransport<Phys>::EnforceVectorPotentialBoundary(IdefixArray4D<re
 
       this->ExchangeAll(Ax1, Ax2, Ax3);
     #endif
+    EnforceEMFBoundaryPeriodic(Ax1, Ax2, Ax3);
   #endif
 
   idfx::popRegion();