Fix BXs consistency accross mpi decomposition

.gitignore

@@ -22,6 +22,9 @@ CMakeFiles
 CMakeCache.txt
 cmake_install.cmake
 Makefile.local
+git-state.txt
+Kokkos_Version_info.*
+
 
 # test artifacts
 test/**/*.o
@@ -36,6 +39,8 @@ test/**/KokkosCore*
 test/**/*.csv
 test/**/*.pyc
 test/**/*.dat
+test/**/cmake_packages*
+test/**/
 
 # machine specific cache and hidden files
 .*

CMakeLists.txt

@@ -99,10 +99,7 @@ if(Idefix_MPI)
   add_compile_definitions("WITH_MPI")
   find_package(MPI REQUIRED)
   target_link_libraries(idefix MPI::MPI_CXX)
-  target_sources(idefix
-    PUBLIC src/mpi.cpp
-    PUBLIC src/mpi.hpp
-  )
+  add_subdirectory(src/mpi)
 endif()
 
 if(Idefix_HDF5)
@@ -234,6 +231,7 @@ target_include_directories(idefix PUBLIC
                            src/gravity
                            src/utils
                            src/utils/iterativesolver
+                           src/mpi
                            src
                            )
 

src/dataBlock/dataBlock.cpp

@@ -367,6 +367,7 @@ real DataBlock::ComputeTimestep() {
 void DataBlock::DeriveVectorPotential() {
   if constexpr(DefaultPhysics::mhd) {
     #ifdef EVOLVE_VECTOR_POTENTIAL
+      hydro->emf->EnforceVectorPotentialBoundary(hydro->Ve);
       hydro->emf->ComputeMagFieldFromA(hydro->Ve, hydro->Vs);
     #endif
   }

src/dataBlock/fargo.cpp

@@ -138,10 +138,15 @@ Fargo::Fargo(Input &input, int nmax, DataBlock *data) {
       for(int i=0 ; i < nvar ; i++) {
         vars.push_back(i);
       }
+      #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
+        const int dirShift = JDIR;
+      #elif GEOMETRY == SPHERICAL
+        const int dirShift = KDIR;
+      #endif
       #if MHD == YES
-        this->mpi.Init(data->mygrid, vars, this->nghost.data(), data->np_int.data(), true);
+        this->mpiExchanger.Init(data->mygrid, dirShift, vars, this->nghost, data->np_int, true);
       #else
-        this->mpi.Init(data->mygrid, vars, this->nghost.data(), data->np_int.data());
+        this->mpiExchanger.Init(data->mygrid, dirShift, vars, this->nghost, data->np_int, false);
       #endif
     }
   #endif

src/dataBlock/fargo.hpp

@@ -11,7 +11,7 @@
 #include <vector>
 #include "idefix.hpp"
 #ifdef WITH_MPI
-  #include "mpi.hpp"
+  #include "exchanger.hpp"
 #endif
 
 #include "physics.hpp"
@@ -63,7 +63,7 @@ class Fargo {
   IdefixArray4D<real> scrhVs;
 
 #ifdef WITH_MPI
-  Mpi mpi;                      // Fargo-specific MPI layer
+  Exchanger mpiExchanger;                      // Fargo-specific MPI layer
 #endif
 
   std::array<int,3> beg;
@@ -290,11 +290,7 @@ void Fargo::StoreToScratch(Fluid<Phys>* hydro) {
   }
   #if WITH_MPI
     if(haveDomainDecomposition) {
-      #if GEOMETRY == CARTESIAN || GEOMETRY == POLAR
-        this->mpi.ExchangeX2(scrhUc, scrhVs);
-      #elif GEOMETRY == SPHERICAL
-        this->mpi.ExchangeX3(scrhUc, scrhVs);
-      #endif
+      this->mpiExchanger.Exchange(scrhUc, scrhVs);
     }
   #endif
 }

src/fluid/boundary/axis.cpp

@@ -22,10 +22,9 @@ void Axis::ShowConfig() {
 }
 
 
-void Axis::SymmetrizeEx1Side(int jref) {
+void Axis::SymmetrizeEx1Side(int jref, IdefixArray3D<real> Ex1) {
 #if DIMENSIONS == 3
 
-  IdefixArray3D<real> Ex1 = this->ex;
   IdefixArray1D<real> Ex1Avg = this->Ex1Avg;
 
   idefix_for("Ex1_ini",0,data->np_tot[IDIR],
@@ -63,9 +62,7 @@ void Axis::SymmetrizeEx1Side(int jref) {
 // Hence, we enforce a regularisation of Ex3 for consistancy.
 
 
-void Axis::RegularizeEx3side(int jref) {
-  IdefixArray3D<real> Ex3 = this->ez;
-
+void Axis::RegularizeEx3side(int jref, IdefixArray3D<real> Ex3) {
   idefix_for("Ex3_Regularise",0,data->np_tot[KDIR],0,data->np_tot[IDIR],
     KOKKOS_LAMBDA(int k,int i) {
       Ex3(k,jref,i) = 0.0;
@@ -133,18 +130,20 @@ void Axis::RegularizeCurrentSide(int side) {
 
 // Average the Emf component along the axis
 
-void Axis::RegularizeEMFs() {
+void Axis::RegularizeEMFs(IdefixArray3D<real> ex,
+                          IdefixArray3D<real> ey,
+                          IdefixArray3D<real> ez) {
   idfx::pushRegion("Axis::RegularizeEMFs");
 
   if(this->axisLeft) {
     int jref = data->beg[JDIR];
-    SymmetrizeEx1Side(jref);
-    RegularizeEx3side(jref);
+    SymmetrizeEx1Side(jref, ex);
+    RegularizeEx3side(jref, ez);
   }
   if(this->axisRight) {
     int jref = data->end[JDIR];
-    SymmetrizeEx1Side(jref);
-    RegularizeEx3side(jref);
+    SymmetrizeEx1Side(jref, ex);
+    RegularizeEx3side(jref, ez);
   }
 
   idfx::popRegion();

src/fluid/boundary/axis.hpp

@@ -26,15 +26,16 @@ class Axis {
  public:
   template <typename Phys>
   explicit Axis(Boundary<Phys> *);  // Initialisation
-  void RegularizeEMFs();                 // Regularize the EMF sitting on the axis
+  void RegularizeEMFs(IdefixArray3D<real>, IdefixArray3D<real>, IdefixArray3D<real>);
+                      // 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 RegularizeBX2s();               // Regularize BX2s on the axis
   void ShowConfig();
 
-
-  void SymmetrizeEx1Side(int);         // Symmetrize on a specific side (internal method)
-  void RegularizeEx3side(int);         // Regularize Ex3 along the axis (internal method)
+  // Internal methods
+  void SymmetrizeEx1Side(int, IdefixArray3D<real>); // Symmetrize on a specific side
+  void RegularizeEx3side(int side, IdefixArray3D<real> ex3); // Regularize Ex3 along the axis
   void RegularizeCurrentSide(int);      // Regularize J along the axis (internal method)
   void FixBx2sAxis(int side);           // Fix BX2s on the axis using the field around it (internal)
   void FixBx2sAxisGhostAverage(int side); //Fix BX2s on the axis using the average of neighbouring
@@ -76,9 +77,6 @@ class Axis {
   IdefixArray1D<int> symmetryVc;
   IdefixArray1D<int> symmetryVs;
 
-  IdefixArray3D<real> ex;
-  IdefixArray3D<real> ey;
-  IdefixArray3D<real> ez;
   IdefixArray4D<real> J;
 
   IdefixArray4D<real> Vc;
@@ -94,11 +92,6 @@ Axis::Axis(Boundary<Phys> *boundary) {
   Vc = boundary->Vc;
   Vs = boundary->Vs;
   J = boundary->fluid->J;
-  if constexpr(Phys::mhd) {
-    ex = boundary->fluid->emf->ex;
-    ey = boundary->fluid->emf->ey;
-    ez = boundary->fluid->emf->ez;
-  }
 
   data = boundary->data;
   haveMHD = Phys::mhd;