diff --git a/src/fluid/boundary/axis.cpp b/src/fluid/boundary/axis.cpp
index 3ad67f0a..e447aefc 100644
--- a/src/fluid/boundary/axis.cpp
+++ b/src/fluid/boundary/axis.cpp
@@ -358,65 +358,28 @@ void Axis::EnforceAxisBoundary(int side) {
 
 // Reconstruct Bx2s taking care of the sides where an axis is lying
 
-void Axis::ReconstructBx2s() {
-  idfx::pushRegion("Axis::ReconstructBx2s");
+void Axis::RegularizeBX2s() {
+  idfx::pushRegion("Axis::RegularizeBX2s");
 #if DIMENSIONS >= 2 && MHD == YES
   IdefixArray4D<real> Vs = this->Vs;
-  IdefixArray3D<real> Ax1=data->A[IDIR];
-  IdefixArray3D<real> Ax2=data->A[JDIR];
-  IdefixArray3D<real> Ax3=data->A[KDIR];
-  int nstart = data->beg[JDIR]-1;
-  int nend = data->end[JDIR];
-  int ntot = data->np_tot[JDIR];
-
-  bool haveleft = axisLeft;
-  bool haveright = axisRight;
-
-  if(haveleft) {
-    // This loop is a copy of ReconstructNormalField, with the proper sign when we cross the axis
-    idefix_for("Axis::ReconstructBX2sLeft",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
-      KOKKOS_LAMBDA (int k, int i) {
-        for(int j = nstart ; j>=0 ; j-- ) {
-          Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i)
-                      +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i)  ,
-                                                                                                ,
-                            - ( Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i)
-                                                            - Ax3(k,j,i) * Vs(BX3s,k,j,i) ))));
-        }
-      }
-    );
-  }
-  if(haveright) {
-    // This loop is a copy of ReconstructNormalField, with the proper sign when we cross the axis
-    idefix_for("Axis::ReconstructBX2sRight",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
-      KOKKOS_LAMBDA (int k, int i) {
-        for(int j = nend ; j<ntot ; j++ ) {
-          Vs(BX2s,k,j+1,i) = 1.0 / Ax2(k,j+1,i) * ( Ax2(k,j,i)*Vs(BX2s,k,j,i)
-                      -(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i)  ,
-                                                                                                ,
-                            -  ( Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i)
-                                                            - Ax3(k,j,i) * Vs(BX3s,k,j,i) ))));
-        }
-      }
-    );
-  }
 
     // Set BX2s on the axis to the average of the two agacent cells
     // This is required since Bx2s on the axis is not evolved since
     // there is no circulation around it
 
-
     int jright = data->end[JDIR];
     int jleft = data->beg[JDIR];
     if(isTwoPi) {
       #ifdef AXIS_BX2S_USE_ATHENA_REGULARISATION
-        if(haveleft) FixBx2sAxisGhostAverage(left);
-        if(haveright) FixBx2sAxisGhostAverage(right);
+        if(axisLeft) FixBx2sAxisGhostAverage(left);
+        if(axisRight) FixBx2sAxisGhostAverage(right);
       #else
-        if(haveleft) FixBx2sAxis(left);
-        if(haveright) FixBx2sAxis(right);
+        if(axisLeft) FixBx2sAxis(left);
+        if(axisRight) FixBx2sAxis(right);
       #endif
     } else {
+      bool haveleft = axisLeft;
+      bool haveright = axisRight;
       idefix_for("Axis:BoundaryAvg",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
             KOKKOS_LAMBDA (int k, int i) {
               if(haveleft) {
diff --git a/src/fluid/boundary/axis.hpp b/src/fluid/boundary/axis.hpp
index a1c3054b..b7469786 100644
--- a/src/fluid/boundary/axis.hpp
+++ b/src/fluid/boundary/axis.hpp
@@ -29,7 +29,7 @@ class Axis {
   void RegularizeEMFs();                 // Regularize the EMF sitting on the axis
   void RegularizeCurrent();             // Regularize the currents along the axis
   void EnforceAxisBoundary(int side);   // Enforce the boundary conditions (along X2)
-  void ReconstructBx2s();               // Reconstruct BX2s in the ghost zone using divB=0
+  void RegularizeBX2s();               // Regularize BX2s on the axis
   void ShowConfig();
 
 
@@ -42,6 +42,8 @@ class Axis {
   void ExchangeMPI(int side);           // Function has to be public for GPU, but its technically
                                         // a private function
 
+  bool haveLeftAxis() const { return axisLeft; }  ///< Check if the axis is on the left
+  bool haveRightAxis() const { return axisRight; } ///< Check if the axis is on the right
 
  private:
   bool isTwoPi = false;
@@ -147,7 +149,6 @@ Axis::Axis(Boundary<Phys> *boundary) {
   Kokkos::deep_copy(symmetryVc, symmetryVcHost);
 
   if constexpr(Phys::mhd) {
-    idfx::cout << "Phys MHD" << std::endl;
     symmetryVs = IdefixArray1D<int>("Axis:SymmetryVs",DIMENSIONS);
     IdefixArray1D<int>::HostMirror symmetryVsHost = Kokkos::create_mirror_view(symmetryVs);
     // Init the array
diff --git a/src/fluid/boundary/boundary.hpp b/src/fluid/boundary/boundary.hpp
index 86909936..897c1ce1 100644
--- a/src/fluid/boundary/boundary.hpp
+++ b/src/fluid/boundary/boundary.hpp
@@ -112,6 +112,8 @@ class Boundary {
   IdefixArray4D<real> Vs; ///< reference to face-centered array that we should sync
   std::unique_ptr<Axis> axis; ///< Axis object, initialised if needed.
   bool haveAxis{false};
+  bool haveLeftAxis{false};  ///< True if the left boundary is an axis
+  bool haveRightAxis{false}; ///< True if the right boundary is an axis
 
  private:
   friend class Axis;
@@ -137,6 +139,8 @@ Boundary<Phys>::Boundary(Fluid<Phys>* fluid) {
   if(data->haveAxis) {
     this->axis = std::make_unique<Axis>(this);
     this->haveAxis = true;
+    this->haveLeftAxis = axis->haveLeftAxis();
+    this->haveRightAxis = axis->haveRightAxis();
   }
 
 
@@ -464,30 +468,34 @@ void Boundary<Phys>::ReconstructNormalField(int dir) {
   if(dir==JDIR) {
     nstart = data->beg[JDIR]-1;
     nend = data->end[JDIR];
-    if(!this->haveAxis) {
-      idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
-        KOKKOS_LAMBDA (int k, int i) {
-          if(reconstructLeft) {
-            for(int j = nstart ; j>=0 ; j-- ) {
-              Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i)
-                      +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i)  ,
-                                                                                                ,
-                            + Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) )));
-            }
+    #if DIMENSIONS == 3
+      const int signLeft = haveLeftAxis ? -1 : 1; // left axis is in the negative direction
+      const int signRight = haveRightAxis ? -1 : 1; // right axis is in the negative direction
+    #endif
+
+    idefix_for("ReconstructBX2s",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
+      KOKKOS_LAMBDA (int k, int i) {
+        if(reconstructLeft) {
+          for(int j = nstart ; j>=0 ; j-- ) {
+            Vs(BX2s,k,j,i) = 1.0 / Ax2(k,j,i) * ( Ax2(k,j+1,i)*Vs(BX2s,k,j+1,i)
+                    +(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i)  ,
+                                                                                              ,
+                    + signLeft*(Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) ))));
           }
-          if(reconstructRight) {
-            for(int j = nend ; j<nx2 ; j++ ) {
-              Vs(BX2s,k,j+1,i) = 1.0 / Ax2(k,j+1,i) * ( Ax2(k,j,i)*Vs(BX2s,k,j,i)
-                       -(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i)  ,
-                                                                                                ,
-                            + Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) )));
-            }
+        }
+        if(reconstructRight) {
+          for(int j = nend ; j<nx2 ; j++ ) {
+            Vs(BX2s,k,j+1,i) = 1.0 / Ax2(k,j+1,i) * ( Ax2(k,j,i)*Vs(BX2s,k,j,i)
+                      -(D_EXPAND( Ax1(k,j,i+1) * Vs(BX1s,k,j,i+1) - Ax1(k,j,i) * Vs(BX1s,k,j,i)  ,
+                                                                                              ,
+                    + signRight*(Ax3(k+1,j,i) * Vs(BX3s,k+1,j,i) - Ax3(k,j,i) * Vs(BX3s,k,j,i) ))));
           }
         }
-      );
-    } else {
+      }
+    );
+    if(haveAxis) {
       // We have an axis, ask myAxis to do that job for us
-      axis->ReconstructBx2s();
+      axis->RegularizeBX2s();
     }
   }
 #endif
@@ -686,7 +694,7 @@ void Boundary<Phys>::EnforceReflective(int dir, BoundarySide side ) {
           const int jref = (dir==JDIR) ? 2*(jghost + side*nxj) - j - 1 : j;
           const int kref = (dir==KDIR) ? 2*(kghost + side*nxk) - k - 1 : k;
 
-          const int sign = (n == VX1+dir) ? -1.0 : 1.0;
+          const int sign = (n == VX1+dir || (n >= BX1 && n != BX1+dir)) ? -1.0 : 1.0;
 
           Vc(n,k,j,i) = sign * Vc(n,kref,jref,iref);
         });