diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 76059965..fc9ba803 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -21,6 +21,44 @@ set(SOURCES
 )
 
 add_library(polyMPO-core ${SOURCES})
+
+# Default ON only for Cray GPU builds (Kokkos CUDA + Cray MPICH GTL available).
+# Can still be overridden with -DpolyMPO_ENABLE_GPU_AWARE_MPI=ON/OFF.
+set(polyMPO_GPU_AWARE_MPI_DEFAULT OFF)
+if(Kokkos_ENABLE_CUDA
+   AND DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80}
+   AND DEFINED ENV{PE_MPICH_GTL_LIBS_nvidia80})
+  set(polyMPO_GPU_AWARE_MPI_DEFAULT ON)
+endif()
+
+option(polyMPO_ENABLE_GPU_AWARE_MPI
+  "Use GPU-aware (CUDA-aware) MPI for halo exchange"
+  ${polyMPO_GPU_AWARE_MPI_DEFAULT})
+message(STATUS "polyMPO_ENABLE_GPU_AWARE_MPI: ${polyMPO_ENABLE_GPU_AWARE_MPI}")
+
+if(polyMPO_ENABLE_GPU_AWARE_MPI)
+  if(NOT Kokkos_ENABLE_CUDA)
+    message(FATAL_ERROR
+      "polyMPO_ENABLE_GPU_AWARE_MPI=ON requires Kokkos built with CUDA.")
+  endif()
+
+  if(DEFINED ENV{PE_MPICH_GTL_DIR_nvidia80}
+     AND DEFINED ENV{PE_MPICH_GTL_LIBS_nvidia80})
+    target_link_options(polyMPO-core PUBLIC
+      "$ENV{PE_MPICH_GTL_DIR_nvidia80}"
+      "$ENV{PE_MPICH_GTL_LIBS_nvidia80}"
+    )
+    message(STATUS "polyMPO: enabling Cray MPICH CUDA GTL linkage")
+  else()
+    message(FATAL_ERROR
+      "Cray MPICH CUDA GTL environment variables are unavailable. "
+      "Load cudatoolkit and craype-accel-nvidia80, or set "
+      "CRAY_ACCEL_TARGET=nvidia80 before configuring.")
+  endif()
+
+  target_compile_definitions(polyMPO-core PUBLIC GPU_AWARE_MPI)
+endif()
+
 set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD "17")
 set_property(TARGET polyMPO-core PROPERTY CXX_STANDARD_REQUIRED ON)
 set_property(TARGET polyMPO-core PROPERTY CXX_EXTENSIONS OFF)
diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp
index 3801dcd9..86f3cfd1 100644
--- a/src/pmpo_MPMesh.cpp
+++ b/src/pmpo_MPMesh.cpp
@@ -166,12 +166,12 @@ void MPMesh::calculateStressDivergence(){
   };
   p_MPs->parallel_for(stress_div, " stress_div_assembly");
   Kokkos::fence();
-  pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds()); 
+  pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds());
 
   timer.reset();
   if(numProcsTot>1){ 
     //Takes contribution of halo vertices and adds it in owner procs
-    communicate_and_take_halo_contributions1(stress_divUV, numVertices, 2, 0, 0);
+    communicate_and_take_halo_contributions_gpu_aware(stress_divUV, numVertices, 2, 0, 0);
     //Transfer the correct values at owned vertices to halo vertices
     //communicate_and_take_halo_contributions(stress_divUV, numVertices, 2, 1, 1);
   }
diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp
index d3aac1b7..2cf71ede 100644
--- a/src/pmpo_MPMesh.hpp
+++ b/src/pmpo_MPMesh.hpp
@@ -4,6 +4,14 @@
 #include "pmpo_utils.hpp"
 #include "pmpo_mesh.hpp"
 #include "pmpo_materialPoints.hpp"
+#include <cstdlib>
+#include <iostream>
+#include <utility>
+#include <stdexcept>
+#include <string>
+#ifdef KOKKOS_ENABLE_CUDA
+#include <cuda_runtime.h>
+#endif
 
 namespace polyMPO{
 
@@ -28,7 +36,7 @@ class MPMesh{
     int numOwnersTot, numHalosTot;
     std::vector<int> numOwnersOnOtherProcs;
     std::vector<int> numHalosOnOtherProcs;
-    std::vector<int>haloOwnerProcs;
+    std::vector<int> haloOwnerProcs;
     std::vector<std::vector<int>> haloOwnerLocalIDs;
     std::vector<std::vector<int>> ownerOwnerLocalIDs;
     std::vector<std::vector<int>> ownerHaloLocalIDs;
@@ -39,7 +47,7 @@ class MPMesh{
     
     //Now Kokkos views are made 1D
     template <typename ViewType>
-    void communicate_and_take_halo_contributions1(
+    void communicate_and_take_halo_contributions_staged(
         const ViewType& meshField,
         int nEntities,
         int numEntries,
@@ -58,10 +66,10 @@ class MPMesh{
       std::vector<std::vector<int>>    recvIDVec;
       std::vector<std::vector<double>> recvDataVec;
       pumipic::RecordTime("SD: Recv Vec Allocation-" + std::to_string(self), timer.seconds());
-      
+
       timer.reset();
-      //communicateFields1(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec);
-      communicateFields1(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec);
+      //communicateFieldsFromHostView(fieldData1, nEntities, numEntries, mode, recvIDVec, recvDataVec);
+      communicateFieldsFromHostView(reconVals_host, nEntities, numEntries, mode, recvIDVec, recvDataVec);
       pumipic::RecordTime("SD: IP Comm-" + std::to_string(self), timer.seconds());
 
       timer.reset();
@@ -112,7 +120,7 @@ class MPMesh{
         assert(recvDataVec[i].size() == recvIDVec[i].size() * numEntries);
       }
       pumipic::RecordTime("SD: Copy CPU-GPU2-" + std::to_string(self), timer.seconds());
- 
+
       //Take contributions from other procs
       timer.reset();
       Kokkos::parallel_for("halo contribution", recvIDGPU.size(), KOKKOS_LAMBDA(const int i){
@@ -128,11 +136,10 @@ class MPMesh{
 
 
     void communicateFields(const std::vector<std::vector<double>>& fieldData, const int numEntities, const int numEntries, int mode,
-                              std::vector<std::vector<int>>& recvIDVec,  std::vector<std::vector<double>>& recvDataVec);
-   
+                             std::vector<std::vector<int>>& recvIDVec, std::vector<std::vector<double>>& recvDataVec);
 
     template <class ViewType> 
-    void communicateFields1(
+    void communicateFieldsFromHostView(
         const ViewType& fieldData, 
         const int numEntities, const int numEntries, int mode,
         std::vector<std::vector<int>>& recvIDVec,
@@ -290,9 +297,400 @@ class MPMesh{
     void calculateStrain();
     void calculateStress(const int constitutive_relation);
     void calculateStressDivergence();
+
+    //Prints the selected halo-exchange path once, on rank 0
+    bool haloExchangePathReported = false;
+    void reportHaloExchangePath(const char* path){
+      if(haloExchangePathReported) return;
+      haloExchangePathReported = true;
+      int self;
+      MPI_Comm_rank(p_MPs->getMPIComm(), &self);
+      if(self == 0){
+        std::cout << "polyMPO halo exchange: " << path << std::endl;
+      }
+    }
+
+#ifdef GPU_AWARE_MPI
+
+#ifdef KOKKOS_ENABLE_CUDA
+    //Device buffer allocated with plain cudaMalloc and exposed as an
+    //unmanaged Kokkos::View. Used for buffers passed directly to MPI.
+    //
+    //Kokkos (>= 4.2, CUDA_MALLOC_ASYNC=ON by default) allocates Views from
+    //cudaMallocAsync memory pools, which cuIpcGetMemHandle rejects
+    //(CUDA_ERROR_INVALID_VALUE). Cray MPICH/GTL relies on CUDA IPC for
+    //intra-node GPU-to-GPU transfers, so these buffers bypass the Kokkos
+    //allocator regardless of how Kokkos was built.
+    template <typename T>
+    struct RawCudaMPIBuffer{
+      T* ptr = nullptr;
+      size_t count = 0;
+
+      void allocate(size_t n){
+        free();
+        count = n;
+        if(n > 0){
+          cudaError_t err = cudaMalloc(&ptr, n * sizeof(T));
+          if(err != cudaSuccess){
+            throw std::runtime_error(std::string("RawCudaMPIBuffer: cudaMalloc failed: ") + cudaGetErrorString(err));
+          }
+        }
+      }
+
+      void free(){
+        if(ptr != nullptr){ cudaFree(ptr); ptr = nullptr; }
+        count = 0;
+      }
+
+      T* data() const{ return ptr; }
+      size_t size() const{ return count; }
+      Kokkos::View<T*, Kokkos::CudaSpace, Kokkos::MemoryUnmanaged> view() const{
+        return Kokkos::View<T*, Kokkos::CudaSpace, Kokkos::MemoryUnmanaged>(ptr, count);
+      }
+
+      RawCudaMPIBuffer() = default;
+      ~RawCudaMPIBuffer(){ free(); }
+      RawCudaMPIBuffer(const RawCudaMPIBuffer&) = delete;
+      RawCudaMPIBuffer& operator=(const RawCudaMPIBuffer&) = delete;
+
+      RawCudaMPIBuffer(RawCudaMPIBuffer&& other) noexcept{
+        ptr = other.ptr; count = other.count;
+        other.ptr = nullptr; other.count = 0;
+      }
+
+      RawCudaMPIBuffer& operator=(RawCudaMPIBuffer&& other) noexcept{
+        if(this != &other){
+          free();
+          ptr = other.ptr; count = other.count;
+          other.ptr = nullptr; other.count = 0;
+        }
+        return *this;
+      }
+    };
+
+    using CudaAwareMPIIntBuffer = RawCudaMPIBuffer<int>;
+    using CudaAwareMPIDoubleBuffer = RawCudaMPIBuffer<double>;
+#else
+    //Non-CUDA backend: plain Kokkos::View with the same interface as RawCudaMPIBuffer
+    template <typename T>
+    struct KokkosMPIBuffer{
+      Kokkos::View<T*> v;
+
+      void allocate(size_t n){
+        v = Kokkos::View<T*>("cudaAwareMPIBuffer_batched", n);
+      }
+
+      T* data() const{ return v.data(); }
+      size_t size() const{ return v.extent(0); }
+      Kokkos::View<T*> view() const{ return v; }
+    };
+
+    using CudaAwareMPIIntBuffer = KokkosMPIBuffer<int>;
+    using CudaAwareMPIDoubleBuffer = KokkosMPIBuffer<double>;
+#endif
+
+    bool cudaAwareMPICacheValid = true;
+    bool cudaAwareMPIDisabled = false;
+    bool mpiGpuSupportChecked = false;
+    bool mpiGpuSupportEnabled = false;
+
+    struct CudaAwareMPIFieldCache{
+      bool valid = false;
+      int cachedNumProcs = -1;
+
+      std::vector<int> sendCounts;
+      std::vector<int> recvCounts;
+      std::vector<int> sendOffsets; //prefix sum of sendCounts, in entities
+      std::vector<int> recvOffsets; //prefix sum of recvCounts, in entities
+
+      int totalSendCount = 0;
+      int totalRecvCount = 0;
+
+      //Batched buffers for all neighbors. Per-proc slices are
+      //[offset, offset + count) for IDs and
+      //[offset * numEntries, (offset + count) * numEntries) for data.
+      CudaAwareMPIIntBuffer sendEntityGPU;
+      CudaAwareMPIIntBuffer recvIDGPU;
+      CudaAwareMPIDoubleBuffer sendDataGPU;
+      CudaAwareMPIDoubleBuffer recvDataGPU;
+    };
+
+    std::map<std::pair<int, int>, CudaAwareMPIFieldCache> cudaAwareMPICaches;
+
+    //The halo-exchange path follows MPICH_GPU_SUPPORT_ENABLED:
+    //  1                                  -> GPU-aware MPI
+    //  anything else (0, unset, yes, ...) -> CPU-staged path
+    //MPI must not be given device pointers unless GPU support is enabled;
+    //doing so crashes instead of returning an MPI error.
+    bool mpiGpuSupport(){
+      if(!mpiGpuSupportChecked){
+        const char* value = std::getenv("MPICH_GPU_SUPPORT_ENABLED");
+        mpiGpuSupportEnabled = value != nullptr && std::string(value) == "1";
+        mpiGpuSupportChecked = true;
+      }
+      return mpiGpuSupportEnabled;
+    }
+
+    //GPU-aware MPI path: field data is sent/received directly from GPU
+    //buffers. Receive IDs come from the fixed halo/owner mapping and are
+    //not exchanged. All neighbors share one batched buffer per direction,
+    //packed/unpacked with a single kernel each.
+    //
+    //Metadata and buffers are cached per (mode, numEntries). If the
+    //communication pattern changes, clear cudaAwareMPICaches first.
+    template <typename ViewType>
+    void communicate_and_take_halo_contributions_gpu_aware(
+        const ViewType& meshField,
+        int nEntities,
+        int numEntries,
+        int mode,
+        int op){
+
+      int self, numProcsTot;
+      MPI_Comm comm = p_MPs->getMPIComm();
+      MPI_Comm_rank(comm, &self);
+      MPI_Comm_size(comm, &numProcsTot);
+
+      assert(mode == 0 || mode == 1);
+      assert(op == 0 || op == 1);
+      assert(nEntities == numOwnersTot + numHalosTot);
+
+      if(cudaAwareMPIDisabled || !mpiGpuSupport()){
+        reportHaloExchangePath("CPU-staged (MPICH_GPU_SUPPORT_ENABLED is not 1)");
+        communicate_and_take_halo_contributions_staged(meshField, nEntities, numEntries, mode, op);
+        return;
+      }
+
+      reportHaloExchangePath("GPU-aware MPI");
+
+      if(!cudaAwareMPICacheValid){
+        cudaAwareMPICaches.clear();
+        cudaAwareMPICacheValid = true;
+      }
+
+      auto& cudaAwareCache = cudaAwareMPICaches[std::make_pair(mode, numEntries)];
+      const bool needRebuild = (!cudaAwareCache.valid) || (cudaAwareCache.cachedNumProcs != numProcsTot);
+
+      if(needRebuild){
+        cudaAwareCache.cachedNumProcs = numProcsTot;
+        cudaAwareCache.sendCounts.assign(numProcsTot, 0);
+        cudaAwareCache.recvCounts.assign(numProcsTot, 0);
+        cudaAwareCache.sendOffsets.assign(numProcsTot, 0);
+        cudaAwareCache.recvOffsets.assign(numProcsTot, 0);
+
+        for(int proc = 0; proc < numProcsTot; proc++){
+          if(proc == self) continue;
+          if(mode == 0){
+            cudaAwareCache.sendCounts[proc] = numOwnersOnOtherProcs[proc];
+            cudaAwareCache.recvCounts[proc] = numHalosOnOtherProcs[proc];
+          }
+          else{
+            cudaAwareCache.sendCounts[proc] = numHalosOnOtherProcs[proc];
+            cudaAwareCache.recvCounts[proc] = numOwnersOnOtherProcs[proc];
+          }
+        }
+
+        int totalSend = 0;
+        int totalRecv = 0;
+        for(int proc = 0; proc < numProcsTot; proc++){
+          cudaAwareCache.sendOffsets[proc] = totalSend;
+          totalSend += cudaAwareCache.sendCounts[proc];
+          cudaAwareCache.recvOffsets[proc] = totalRecv;
+          totalRecv += cudaAwareCache.recvCounts[proc];
+        }
+        cudaAwareCache.totalSendCount = totalSend;
+        cudaAwareCache.totalRecvCount = totalRecv;
+
+        cudaAwareCache.sendEntityGPU.allocate(totalSend);
+        cudaAwareCache.sendDataGPU.allocate(totalSend * numEntries);
+        cudaAwareCache.recvIDGPU.allocate(totalRecv);
+        cudaAwareCache.recvDataGPU.allocate(totalRecv * numEntries);
+
+        //Build flattened send-entity list on host, then copy to device
+        if(totalSend > 0){
+          auto sendEntityCPU = Kokkos::View<int*, Kokkos::HostSpace>("sendEntityCPU_batched", totalSend);
+          if(mode == 0){
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+              if(cudaAwareCache.sendCounts[proc] <= 0) continue;
+              assert(haloOwnerLocalIDs[proc].size() == (size_t)cudaAwareCache.sendCounts[proc]);
+            }
+
+            std::vector<int> cursor(cudaAwareCache.sendOffsets);
+            for(int iEnt = 0; iEnt < numHalosTot; iEnt++){
+              int ownerProc = haloOwnerProcs[iEnt];
+              if(ownerProc == self) continue;
+              sendEntityCPU(cursor[ownerProc]) = numOwnersTot + iEnt;
+              cursor[ownerProc]++;
+            }
+
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+              assert(cursor[proc] == cudaAwareCache.sendOffsets[proc] + cudaAwareCache.sendCounts[proc]);
+            }
+          }
+          else{
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+              int sendCount = cudaAwareCache.sendCounts[proc];
+              if(sendCount <= 0) continue;
+              assert(ownerOwnerLocalIDs[proc].size() == (size_t)sendCount);
+              int base = cudaAwareCache.sendOffsets[proc];
+              for(int i = 0; i < sendCount; i++){
+                sendEntityCPU(base + i) = ownerOwnerLocalIDs[proc][i];
+              }
+            }
+          }
+          Kokkos::deep_copy(cudaAwareCache.sendEntityGPU.view(), sendEntityCPU);
+        }
+
+        //Build flattened recv-ID list on host, then copy to device
+        if(totalRecv > 0){
+          auto recvIDCPU = Kokkos::View<int*, Kokkos::HostSpace>("recvIDCPU_batched", totalRecv);
+          if(mode == 0){
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+              int recvCount = cudaAwareCache.recvCounts[proc];
+              if(recvCount <= 0) continue;
+              assert(ownerOwnerLocalIDs[proc].size() == (size_t)recvCount);
+              int base = cudaAwareCache.recvOffsets[proc];
+              for(int i = 0; i < recvCount; i++){
+                recvIDCPU(base + i) = ownerOwnerLocalIDs[proc][i];
+              }
+            }
+          }
+          else{
+            std::vector<int> cursor(cudaAwareCache.recvOffsets);
+            for(int iEnt = 0; iEnt < numHalosTot; iEnt++){
+              int ownerProc = haloOwnerProcs[iEnt];
+              if(ownerProc == self) continue;
+              recvIDCPU(cursor[ownerProc]) = numOwnersTot + iEnt;
+              cursor[ownerProc]++;
+            }
+
+            for(int proc = 0; proc < numProcsTot; proc++){
+              if(proc == self) continue;
+              assert(cursor[proc] == cudaAwareCache.recvOffsets[proc] + cudaAwareCache.recvCounts[proc]);
+            }
+          }
+          Kokkos::deep_copy(cudaAwareCache.recvIDGPU.view(), recvIDCPU);
+        }
+
+        cudaAwareCache.valid = true;
+      }
+
+      //Pack send buffer
+      if(cudaAwareCache.totalSendCount > 0){
+        auto sendEntityGPU = cudaAwareCache.sendEntityGPU.view();
+        auto sendDataGPU = cudaAwareCache.sendDataGPU.view();
+        Kokkos::parallel_for("pack cached gpu-aware mpi send buffer batched", cudaAwareCache.totalSendCount, KOKKOS_LAMBDA(const int i){
+          int entity = sendEntityGPU(i);
+          for(int k = 0; k < numEntries; k++){
+            sendDataGPU(i * numEntries + k) = meshField(entity, k);
+          }
+        });
+      }
+      //MPI is not stream-aware: packing must finish before Isend
+      Kokkos::fence();
+
+      std::vector<MPI_Request> requests;
+      requests.reserve(2 * numProcsTot);
+      int mpiError = MPI_SUCCESS;
+
+      //Post Irecv/Isend per neighbor
+      for(int proc = 0; proc < numProcsTot; proc++){
+        if(proc == self) continue;
+
+        if(cudaAwareCache.recvCounts[proc] > 0){
+          MPI_Request reqData;
+          double* recvPtr = cudaAwareCache.recvDataGPU.data() + (size_t)cudaAwareCache.recvOffsets[proc] * numEntries;
+          mpiError = MPI_Irecv(recvPtr, cudaAwareCache.recvCounts[proc] * numEntries, MPI_DOUBLE, proc, 2, comm, &reqData);
+          if(mpiError != MPI_SUCCESS) break;
+          requests.push_back(reqData);
+        }
+
+        if(cudaAwareCache.sendCounts[proc] > 0){
+          MPI_Request reqData;
+          double* sendPtr = cudaAwareCache.sendDataGPU.data() + (size_t)cudaAwareCache.sendOffsets[proc] * numEntries;
+          mpiError = MPI_Isend(sendPtr, cudaAwareCache.sendCounts[proc] * numEntries, MPI_DOUBLE, proc, 2, comm, &reqData);
+          if(mpiError != MPI_SUCCESS) break;
+          requests.push_back(reqData);
+        }
+      }
+
+      if(mpiError == MPI_SUCCESS && !requests.empty()){
+        mpiError = MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+      }
+
+      if(mpiError != MPI_SUCCESS){
+        cudaAwareMPIDisabled = true;
+        if(self == 0){
+          std::cout << "[GPU_AWARE_MPI] Batched device-pointer MPI failed." << std::endl;
+        }
+
+        if(requests.empty()){
+          if(self == 0){
+            std::cout << "[GPU_AWARE_MPI] Falling back to CPU-staged communication." << std::endl;
+          }
+          communicate_and_take_halo_contributions_staged(meshField, nEntities, numEntries, mode, op);
+          return;
+        }
+
+        if(self == 0){
+          std::cout << "[GPU_AWARE_MPI] Failure happened after MPI requests were posted. "
+                    << "Set MPICH_GPU_SUPPORT_ENABLED=0 before running to use the CPU-staged path." << std::endl;
+        }
+        MPI_Abort(comm, mpiError);
+        return;
+      }
+
+      //Unpack received contributions
+      if(cudaAwareCache.totalRecvCount > 0){
+        auto recvIDGPU = cudaAwareCache.recvIDGPU.view();
+        auto recvDataGPU = cudaAwareCache.recvDataGPU.view();
+        if(op == 0){
+          Kokkos::parallel_for("halo add cached gpu-aware mpi batched", cudaAwareCache.totalRecvCount, KOKKOS_LAMBDA(const int i){
+            const int vertex = recvIDGPU(i);
+            for(int k = 0; k < numEntries; k++){
+#ifdef POLYMPO_ASSUME_UNIQUE_HALO_CONTRIBS
+              meshField(vertex, k) += recvDataGPU(i * numEntries + k);
+#else
+              Kokkos::atomic_add(&meshField(vertex, k), recvDataGPU(i * numEntries + k));
+#endif
+            }
+          });
+        }
+        else{
+          Kokkos::parallel_for("halo assign cached gpu-aware mpi batched", cudaAwareCache.totalRecvCount, KOKKOS_LAMBDA(const int i){
+            const int vertex = recvIDGPU(i);
+            for(int k = 0; k < numEntries; k++){
+              meshField(vertex, k) = recvDataGPU(i * numEntries + k);
+            }
+          });
+        }
+      }
+      Kokkos::fence();
+    }
+
+#else
+
+    //GPU_AWARE_MPI not defined: use the CPU-staged path
+    template <typename ViewType>
+    void communicate_and_take_halo_contributions_gpu_aware(
+        const ViewType& meshField,
+        int nEntities,
+        int numEntries,
+        int mode,
+        int op){
+
+      reportHaloExchangePath("CPU-staged (built without GPU_AWARE_MPI)");
+      communicate_and_take_halo_contributions_staged(meshField, nEntities, numEntries, mode, op);
+    }
+
+#endif
+
 };
 
 }//namespace polyMPO end
 
 #endif
-
diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp
index 9f184e05..c26d3899 100644
--- a/src/pmpo_MPMesh_assembly.hpp
+++ b/src/pmpo_MPMesh_assembly.hpp
@@ -161,14 +161,14 @@ void MPMesh::reconstruct_coeff_full(){
   int mode = 0;
   int op = 0;
   if (numProcsTot >1){
-    communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
     mode=1; 
     op=1;
-    communicate_and_take_halo_contributions1(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
+    communicate_and_take_halo_contributions_gpu_aware(vtxMatrices, numVertices, numEntriesMatrix, mode, op);
   }
   pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds());
  
-  //Stroe the 1st matrix element
+  //Store the 1st matrix element
   Kokkos::View<double*>vtxMatrixMass_l("vtxMass", numVertices);
   Kokkos::parallel_for("storeMatrixMass", numVertices, KOKKOS_LAMBDA(const int vtx){
     vtxMatrixMass_l(vtx) = vtxMatrices(vtx, 0);
@@ -380,7 +380,7 @@ void MPMesh::assemblyVtx1(){
 
   timer.reset();
   if(numProcsTot>1){ 
-    communicate_and_take_halo_contributions1(meshField, numVertices, numEntries, 0, 0);
+    communicate_and_take_halo_contributions_gpu_aware(meshField, numVertices, numEntries, 0, 0);
   }
   pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds());
 }
diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index 53c330e1..2ee56258 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -1605,11 +1605,19 @@ void polympo_set_free_slip_bc_f(MPMesh_ptr p_mpmesh){
 }
 
 void polympo_set_halo_vel_from_owner_f(MPMesh_ptr p_mpmesh){
+
+  int numProcsTot;
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  MPI_Comm comm = p_MPs->getMPIComm();
+  MPI_Comm_size(comm, &numProcsTot);
+  if(numProcsTot == 1) return;
+
   auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
   int numVertices = p_mesh->getNumVertices();
   auto vtxFieldVel = p_mesh->getMeshField<polyMPO::MeshF_Vel>();
-  mpMesh->communicate_and_take_halo_contributions1(vtxFieldVel, numVertices, 2, 1, 1); 
+
+  mpMesh->communicate_and_take_halo_contributions_gpu_aware(vtxFieldVel, numVertices, 2, 1, 1);
 }
 
 //Advection Calcualtions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 3183cdc5..346c01b4 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -61,6 +61,7 @@ pmpo_add_exe(testReconstruction testReconstruction.cpp)
 pmpo_add_exe(testTracking testTracking.cpp)
 pmpo_add_exe(testRebuild testRebuild.cpp)
 pmpo_add_exe(timeAssmblyWachspress testTiming.cpp)
+pmpo_add_exe(testGPUAwareHaloExchange testGPUAwareHaloExchange.cpp)
 pmpo_add_fortran_exe(testFortranInit testFortranInit.f90)
 pmpo_add_fortran_exe(testFortran testFortran.f90)
 pmpo_add_fortran_exe(testFortranReadMPAS testFortranReadMPAS.f90)
@@ -91,6 +92,8 @@ pmpo_add_test(testFortranInit 1 ./testFortranInit)
 pmpo_add_test(testFortran 1 ./testFortran)
 pmpo_add_test(testFortranMPMeshModule 1 ./testFortranMPMeshModule)
 pmpo_add_test(testFortranMPAppIDs 1 ./testFortranMPAppIDs)
+pmpo_add_test(test_gpu_aware_halo_exchange 4 ./testGPUAwareHaloExchange)
+set_tests_properties(test_gpu_aware_halo_exchange PROPERTIES SKIP_RETURN_CODE 77)
 
 #set NC file for test
 set(TEST_NC_FILE_PLANAR "${CMAKE_SOURCE_DIR}/test/sample_mpas_meshes/planar_nonuniform_cvt_for_square_673elms.nc")
diff --git a/test/testGPUAwareHaloExchange.cpp b/test/testGPUAwareHaloExchange.cpp
new file mode 100644
index 00000000..127fa303
--- /dev/null
+++ b/test/testGPUAwareHaloExchange.cpp
@@ -0,0 +1,286 @@
+#include "pmpo_MPMesh.hpp"
+#include "pmpo_createTestMPMesh.hpp"
+
+#include <mpi.h>
+#include <Kokkos_Core.hpp>
+
+#include <iostream>
+
+int main(int argc, char** argv)
+{
+    MPI_Init(&argc, &argv);
+    Kokkos::initialize(argc, argv);
+
+    int testResult = 0;
+
+    {
+        int rank = -1;
+        int size = -1;
+
+        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+        MPI_Comm_size(MPI_COMM_WORLD, &size);
+
+        if (rank == 0) {
+            std::cout
+                << "GPU-aware halo exchange test running with "
+                << size << " MPI ranks."
+                << std::endl;
+        }
+
+        // This test uses the following communication pattern:
+        //
+        // Rank 0 ----\
+        // Rank 1 ----- > Rank 3
+        // Rank 2 ----/
+        //                |
+        //                v
+        //              Rank 0
+        //
+        // Therefore:
+        // Rank 0 has one halo owned by Rank 3.
+        // Rank 1 has no halos.
+        // Rank 2 has no halos.
+        // Rank 3 has three halos owned by Ranks 0, 1, and 2.
+
+        if (size != 4) {
+            if (rank == 0) {
+                std::cerr
+                    << "This test requires exactly 4 MPI ranks (got "
+                    << size << "); skipping."
+                    << std::endl;
+            }
+            Kokkos::finalize();
+            MPI_Finalize();
+            return 77;
+        }
+
+        // Create the existing polyMPO test mesh and MPMesh.
+
+        polyMPO::Mesh* mesh = polyMPO::initTestMesh(1, 1);
+
+        polyMPO::MPMesh mpMesh = polyMPO::initTestMPMesh(mesh, 1);
+
+        // initTestMPMesh() creates MaterialPoints, but the test helper
+        // does not initialize its MPI communicator.
+        mpMesh.p_MPs->setMPIComm(MPI_COMM_WORLD);
+
+        const int nVertices = mesh->getNumVertices();
+
+        // startCommunication() expects local vertices to be ordered:
+
+        polyMPO::IntView owningProcVertex("testOwningProcVertex", nVertices);
+
+        polyMPO::IntView globalVtx("testGlobalVtx",nVertices);
+
+        auto owningProcHost = Kokkos::create_mirror_view(owningProcVertex);
+
+        auto globalVtxHost = Kokkos::create_mirror_view(globalVtx);
+
+        // Number of locally-owned vertices differs by rank.
+        // Rank 0: 18 owners + 1 halo
+        // Rank 1: 19 owners + 0 halos
+        // Rank 2: 19 owners + 0 halos
+        // Rank 3: 16 owners + 3 halos
+
+        int numLocalOwners = nVertices;
+
+        if (rank == 0) {
+            numLocalOwners = nVertices - 1;
+        }
+        else if (rank == 3) {
+            numLocalOwners = nVertices - 3;
+        }
+
+        // Assign locally-owned vertices.
+        // Rank 0 global IDs:   0,   1,   2, ...
+        // Rank 1 global IDs: 100, 101, 102, ...
+        // Rank 2 global IDs: 200, 201, 202, ...
+        // Rank 3 global IDs: 300, 301, 302, ...
+
+        for (int i = 0; i < numLocalOwners; ++i) {
+            owningProcHost(i) = rank;
+            globalVtxHost(i) = rank * 100 + i;
+        }
+
+        // Define halo vertices.
+
+        if (rank == 0) {
+
+            // Rank 0 receives Rank 3 owner vertex 0.
+            // Rank 3 owner vertex 0 has global ID 300.
+
+            owningProcHost(nVertices - 1) = 3;
+            globalVtxHost(nVertices - 1) = 300;
+        }
+        else if (rank == 3) {
+
+            // Rank 3 receives Rank 0 owner vertex 0.
+            owningProcHost(nVertices - 3) = 0;
+            globalVtxHost(nVertices - 3) = 0;
+
+            // Rank 3 receives Rank 1 owner vertex 0.
+            owningProcHost(nVertices - 2) = 1;
+            globalVtxHost(nVertices - 2) = 100;
+
+            // Rank 3 receives Rank 2 owner vertex 0.
+            owningProcHost(nVertices - 1) = 2;
+            globalVtxHost(nVertices - 1) = 200;
+        }
+
+        Kokkos::deep_copy(owningProcVertex, owningProcHost);
+
+        Kokkos::deep_copy(globalVtx, globalVtxHost);
+
+        mesh->setOwningProcVertex(owningProcVertex);
+        mesh->setVtxGlobal(globalVtx);
+
+        // Build owner/halo communication metadata
+        mpMesh.startCommunication();
+
+        std::cout
+            << "Rank " << rank
+            << ": owners = " << mpMesh.numOwnersTot
+            << ", halos = " << mpMesh.numHalosTot
+            << std::endl;
+
+        // Create one scalar field value per vertex.
+
+        const int nEntities = mpMesh.numOwnersTot + mpMesh.numHalosTot;
+
+        const int numEntries = 1;
+
+        Kokkos::View<double**> field("gpuAwareHaloTestField", nEntities, numEntries);
+
+        const int numOwners = mpMesh.numOwnersTot;
+
+        // Owner values are deterministic:
+        // Rank 0 owner 0 =    0
+        // Rank 1 owner 0 = 1000
+        // Rank 2 owner 0 = 2000
+        // Rank 3 owner 0 = 3000
+        // All halos start at -1.
+
+        Kokkos::parallel_for("initialize_gpu_aware_halo_test", nEntities,KOKKOS_LAMBDA(const int i)
+            {
+                if (i < numOwners) {
+                    field(i, 0) =
+                        1000.0 * static_cast<double>(rank)
+                        + static_cast<double>(i);
+                }
+                else {
+                    field(i, 0) = -1.0;
+                }
+            });
+
+        Kokkos::fence();
+
+        // Perform owner -> halo assignment.
+        // mode = 1 : owner -> halo
+        // op   = 1 : assignment
+
+        const int mode = 1;
+        const int op = 1;
+
+        mpMesh.communicate_and_take_halo_contributions_gpu_aware(field, nEntities, numEntries, mode, op);
+
+        Kokkos::fence();
+
+        // Copy results to host and verify exact halo values.
+
+        auto fieldHost = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), field);
+
+        int localFailures = 0;
+
+        if (rank == 0) {
+
+            // Rank 0 receives Rank 3 owner vertex 0.
+            const int haloIndex = nVertices - 1;
+            const double expectedValue = 3000.0;
+
+            if (fieldHost(haloIndex, 0) != expectedValue) {
+
+                ++localFailures;
+
+                std::cerr
+                    << "Rank 0: halo value = "
+                    << fieldHost(haloIndex, 0)
+                    << ", expected = "
+                    << expectedValue
+                    << std::endl;
+            }
+        }
+        else if (rank == 3) {
+
+            // Rank 3 receives owner vertex 0 from Ranks 0, 1, and 2.
+
+            const int haloFromRank0 = nVertices - 3;
+            const int haloFromRank1 = nVertices - 2;
+            const int haloFromRank2 = nVertices - 1;
+
+            if (fieldHost(haloFromRank0, 0) != 0.0) {
+
+                ++localFailures;
+
+                std::cerr
+                    << "Rank 3: halo from Rank 0 = "
+                    << fieldHost(haloFromRank0, 0)
+                    << ", expected = 0"
+                    << std::endl;
+            }
+
+            if (fieldHost(haloFromRank1, 0) != 1000.0) {
+
+                ++localFailures;
+
+                std::cerr
+                    << "Rank 3: halo from Rank 1 = "
+                    << fieldHost(haloFromRank1, 0)
+                    << ", expected = 1000"
+                    << std::endl;
+            }
+
+            if (fieldHost(haloFromRank2, 0) != 2000.0) {
+
+                ++localFailures;
+
+                std::cerr
+                    << "Rank 3: halo from Rank 2 = "
+                    << fieldHost(haloFromRank2, 0)
+                    << ", expected = 2000"
+                    << std::endl;
+            }
+        }
+
+        int globalFailures = 0;
+
+        MPI_Allreduce(&localFailures, &globalFailures, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+        if (rank == 0) {
+
+            if (globalFailures == 0) {
+                std::cout
+                    << "GPU-aware halo exchange test PASSED."
+                    << std::endl;
+            }
+            else {
+                std::cerr
+                    << "GPU-aware halo exchange test FAILED with "
+                    << globalFailures
+                    << " halo errors."
+                    << std::endl;
+            }
+        }
+
+        if (globalFailures != 0) {
+            testResult = 1;
+        }
+
+        // Do not delete mesh here.
+        // mpMesh owns the Mesh and MaterialPoints objects.
+    }
+
+    Kokkos::finalize();
+    MPI_Finalize();
+
+    return testResult;
+}