feat: add halo exchange into dynamics

Author: TomMelt

dynamics/src/CGDynamicsKernel.cpp

@@ -19,6 +19,10 @@
 
 #include <limits>
 
+#ifdef USE_MPI
+#include "include/Halo.hpp"
+#endif
+
 namespace Nextsim {
 
 template <int DGadvection>
@@ -272,6 +276,14 @@ template <int DGadvection> void CGDynamicsKernel<DGadvection>::prepareIteration(
     Interpolations::DG2CG(*smesh, cgA, data.at(ciceName));
     VectorManipulations::CGAveragePeriodic(*smesh, cgA);
 
+#ifdef USE_MPI
+    // Halo object only depends on shape of the array, so we can share it for all DGVectors of the
+    // same shape
+    Halo halo(cgH);
+    halo.exchangeHalos(cgH);
+    halo.exchangeHalos(cgA);
+#endif
+
     // Reinit the gradient of the sea surface height. Not done by
     // DataMap as seaSurfaceHeight is always dG(0)
     computeGradientOfSeaSurfaceHeight(DynamicsKernel<DGadvection, DGstressComp>::seaSurfaceHeight);

dynamics/src/include/BrittleCGDynamicsKernel.hpp

@@ -15,6 +15,10 @@
 #include "include/constants.hpp"
 #include <cmath>
 
+#ifdef USE_MPI
+#include "include/Halo.hpp"
+#endif
+
 namespace Nextsim {
 
 // The brittle momentum solver for CG velocity fields
@@ -116,6 +120,20 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
     {
         advectDynamicsFields(tst.step.seconds());
 
+        // halo exchange
+        //  - damage, hice, cice
+        //  - s11, s12, s22
+        //  - e11, e12, e22
+        //  - dStressY, dStressX
+        //  - u, v
+        // Note: only need to create one halo object per array type, dimensionality and size.
+#ifdef USE_MPI
+        Halo halo(hice);
+        halo.exchangeHalos(static_cast<DGVector<DGadvection>&>(hice));
+        halo.exchangeHalos(static_cast<DGVector<DGadvection>&>(cice));
+        halo.exchangeHalos(static_cast<DGVector<DGadvection>&>(damage));
+#endif
+
         prepareIteration({ { hiceName, hice }, { ciceName, cice } });
 
         // The timestep for the brittle solver is the solver subtimestep
@@ -128,18 +146,42 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
 
             projectVelocityToStrain();
 
+#ifdef USE_MPI
+            Halo haloDGVector(e11);
+            haloDGVector.exchangeHalos(e11);
+            haloDGVector.exchangeHalos(e12);
+            haloDGVector.exchangeHalos(e22);
+#endif
+
             std::array<std::reference_wrapper<DGVector<DGstressComp>>, N_TENSOR_ELEMENTS> stress
                 = { s11, s12, s22 }; // Call the step function on the StressUpdateStep class
             // Call the step function on the StressUpdateStep class
             stressStep.stressUpdateHighOrder(
                 params, *smesh, stress, { e11, e12, e22 }, hice, cice, deltaT);
 
+#ifdef USE_MPI
+            haloDGVector.exchangeHalos(s11);
+            haloDGVector.exchangeHalos(s12);
+            haloDGVector.exchangeHalos(s22);
+#endif
+
             stressDivergence(); // Compute divergence of stress tensor
 
+#ifdef USE_MPI
+            Halo haloCGVector(dStressX);
+            haloCGVector.exchangeHalos(dStressX);
+            haloCGVector.exchangeHalos(dStressY);
+#endif
+
             updateMomentum(tst);
 
             applyBoundaries();
 
+#ifdef USE_MPI
+            haloCGVector.exchangeHalos(u);
+            haloCGVector.exchangeHalos(v);
+#endif
+
             // Land mask
         }