diff --git a/.github/workflows/clang.yml b/.github/workflows/clang.yml
index c996db63c1..d02342899b 100644
--- a/.github/workflows/clang.yml
+++ b/.github/workflows/clang.yml
@@ -126,6 +126,53 @@ jobs:
         ccache -s
         du -hs ~/.cache/ccache
 
+  tests_cxx20:
+    name: Clang C++20 [tests]
+    runs-on: ubuntu-22.04
+    steps:
+    - uses: actions/checkout@v4
+    - name: Dependencies
+      run: |
+        .github/workflows/dependencies/dependencies.sh
+        .github/workflows/dependencies/dependencies_clang.sh 15
+        .github/workflows/dependencies/dependencies_ccache.sh
+    - name: Set Up Cache
+      uses: actions/cache@v4
+      with:
+        path: ~/.cache/ccache
+        key: ccache-${{ github.workflow }}-${{ github.job }}-git-${{ github.sha }}
+        restore-keys: |
+             ccache-${{ github.workflow }}-${{ github.job }}-git-
+    - name: Build & Install
+      env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wnon-virtual-dtor"}
+      run: |
+        export CCACHE_COMPRESS=1
+        export CCACHE_COMPRESSLEVEL=10
+        export CCACHE_MAXSIZE=100M
+        export CCACHE_EXTRAFILES=${{ github.workspace }}/.clang-tidy
+        export CCACHE_LOGFILE=${{ github.workspace }}/ccache.log.txt
+        ccache -z
+
+        mkdir build
+        cd build
+        cmake ..                                      \
+            -DCMAKE_BUILD_TYPE=RelWithDebInfo         \
+            -DCMAKE_VERBOSE_MAKEFILE=ON               \
+            -DAMReX_EB=ON                             \
+            -DAMReX_ENABLE_TESTS=ON                   \
+            -DAMReX_FORTRAN=OFF                       \
+            -DAMReX_MPI=ON                            \
+            -DCMAKE_CXX_STANDARD=20                   \
+            -DCMAKE_C_COMPILER=$(which clang-15)      \
+            -DCMAKE_CXX_COMPILER=$(which clang++-15)  \
+            -DCMAKE_CXX_COMPILER_LAUNCHER=ccache
+        make -j 4
+
+        ctest --output-on-failure
+
+        ccache -s
+        du -hs ~/.cache/ccache
+
   # Build 2D libamrex with configure
   configure-2d:
     name: Clang NOMPI Release [configure 2D]
diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml
index 4090f9c2ec..e0e7bf6911 100644
--- a/.github/workflows/windows.yml
+++ b/.github/workflows/windows.yml
@@ -87,16 +87,17 @@ jobs:
 
         cmake -S . -B build   `
               -DCMAKE_VERBOSE_MAKEFILE=ON   `
+              -DCMAKE_BUILD_TYPE=RelWithDebInfo `
               -DCMAKE_CXX_FLAGS=" /D_DISABLE_CONSTEXPR_MUTEX_CONSTRUCTOR" `
               -DAMReX_EB=ON                 `
               -DAMReX_ENABLE_TESTS=ON       `
               -DAMReX_FORTRAN=OFF           `
               -DAMReX_MPI=OFF
               #-DCMAKE_CXX_COMPILER_LAUNCHER=ccache
-        cmake --build build --config RelWithDebInfo -j 4
+        cmake --build build -j 4
 
-        cmake --build build --config RelWithDebInfo --target install
-        cmake --build build --config RelWithDebInfo --target test_install
+        cmake --build build --target install
+        cmake --build build --target test_install
 
         #ccache -s
 
diff --git a/Src/AmrCore/AMReX_AmrMesh.H b/Src/AmrCore/AMReX_AmrMesh.H
index f5d49f5c5d..27ba1317ca 100644
--- a/Src/AmrCore/AMReX_AmrMesh.H
+++ b/Src/AmrCore/AMReX_AmrMesh.H
@@ -178,7 +178,7 @@ public:
     void SetGeometry (int lev, const Geometry& geom_in) noexcept;
 
     //! Given domain box, return AMR level.  Return -1 if there is no match.
-    int GetLevel (Box const& domain) noexcept;
+    [[nodiscard]] int GetLevel (Box const& domain) const noexcept;
 
     void ClearDistributionMap (int lev) noexcept;
     void ClearBoxArray (int lev) noexcept;
@@ -201,7 +201,7 @@ public:
     //! Return the largest allowable grid.
     [[nodiscard]] const IntVect& maxGridSize (int lev) const noexcept { return max_grid_size[lev]; }
 
-    [[nodiscard]] bool LevelDefined (int lev) noexcept;
+    [[nodiscard]] bool LevelDefined (int lev) const noexcept;
 
     //! Should we keep the coarser grids fixed (and not regrid those levels) at all?
     [[nodiscard]] bool useFixedCoarseGrids () const noexcept { return use_fixed_coarse_grids; }
@@ -255,7 +255,7 @@ public:
 
     [[nodiscard]] virtual BoxArray GetAreaNotToTag (int /*lev*/) { return BoxArray(); }
 
-    [[nodiscard]] long CountCells (int lev) noexcept;
+    [[nodiscard]] Long CountCells (int lev) const noexcept;
 
     [[nodiscard]] virtual DistributionMapping MakeDistributionMap (int lev, BoxArray const& ba);
 
diff --git a/Src/AmrCore/AMReX_AmrMesh.cpp b/Src/AmrCore/AMReX_AmrMesh.cpp
index 0ed59002f2..57b46ca686 100644
--- a/Src/AmrCore/AMReX_AmrMesh.cpp
+++ b/Src/AmrCore/AMReX_AmrMesh.cpp
@@ -421,9 +421,9 @@ AmrMesh::SetGeometry (int lev, const Geometry& geom_in) noexcept
 }
 
 int
-AmrMesh::GetLevel (Box const& domain) noexcept
+AmrMesh::GetLevel (Box const& domain) const noexcept
 {
-    Box ccdomain = amrex::enclosedCells(domain);
+    Box const& ccdomain = amrex::enclosedCells(domain);
     for (int lev = 0; lev < geom.size(); ++lev) {
         if (geom[lev].Domain() == ccdomain) { return lev; }
     }
@@ -443,7 +443,7 @@ AmrMesh::ClearBoxArray (int lev) noexcept
 }
 
 bool
-AmrMesh::LevelDefined (int lev) noexcept
+AmrMesh::LevelDefined (int lev) const noexcept
 {
     return lev <= max_level && !grids[lev].empty() && !dmap[lev].empty();
 }
@@ -1216,8 +1216,8 @@ AmrMesh::checkInput ()
     }
 }
 
-long
-AmrMesh::CountCells (int lev) noexcept
+Long
+AmrMesh::CountCells (int lev) const noexcept
 {
     return grids[lev].numPts();
 }
diff --git a/Src/AmrCore/AMReX_ErrorList.cpp b/Src/AmrCore/AMReX_ErrorList.cpp
index c9d1a7feef..7e72cd956a 100644
--- a/Src/AmrCore/AMReX_ErrorList.cpp
+++ b/Src/AmrCore/AMReX_ErrorList.cpp
@@ -258,9 +258,9 @@ AMRErrorTag::operator() (TagBoxArray&    tba,
             auto const& tagma = tba.arrays();
             if (m_test == BOX)
             {
-                const auto plo = geom.ProbLoArray();
-                const auto dx  = geom.CellSizeArray();
-                const auto tag_rb = m_info.m_realbox;
+                const auto& plo = geom.ProbLoArray();
+                const auto& dx  = geom.CellSizeArray();
+                const auto& tag_rb = m_info.m_realbox;
                 ParallelFor(tba, [=] AMREX_GPU_DEVICE (int bi, int i, int j, int k) noexcept
                 {
                     GpuArray<Real,AMREX_SPACEDIM> pt
@@ -277,7 +277,7 @@ AMRErrorTag::operator() (TagBoxArray&    tba,
                 auto const& datma   = mf->const_arrays();
                 auto threshold = m_value[level];
                 auto const volume_weighting = m_info.m_volume_weighting;
-                auto geomdata = geom.data();
+                auto const& geomdata = geom.data();
                 auto tag_update = tagval;
                 if (m_info.m_derefine) {
                     tag_update = clearval;
diff --git a/Src/Base/AMReX_BaseFwd.H b/Src/Base/AMReX_BaseFwd.H
index 824a6c504e..622689a2e8 100644
--- a/Src/Base/AMReX_BaseFwd.H
+++ b/Src/Base/AMReX_BaseFwd.H
@@ -1,5 +1,6 @@
 #ifndef AMREX_BASE_FWD_H_
 #define AMREX_BASE_FWD_H_
+#include <AMReX_Config.H>
 
 namespace amrex {
 
@@ -25,7 +26,9 @@ class Box;
 template<int dim>
 class IntVectND;
 using IntVect = IntVectND<AMREX_SPACEDIM>;
-class IndexType;
+template<int dim>
+class IndexTypeND;
+using IndexType = IndexTypeND<AMREX_SPACEDIM>;
 struct Dim3;
 struct XDim3;
 class RealBox;
diff --git a/Src/Base/AMReX_GpuRange.H b/Src/Base/AMReX_GpuRange.H
index ecf9a32fd2..08c77595e7 100644
--- a/Src/Base/AMReX_GpuRange.H
+++ b/Src/Base/AMReX_GpuRange.H
@@ -51,7 +51,7 @@ AMREX_FORCE_INLINE Box at (Box const& b, Long offset) noexcept
                                 static_cast<int>(j),
                                 static_cast<int>(k))};
         iv += b.smallEnd();
-        return Box(iv,iv,b.type());
+        return Box(iv,iv,b.ixType());
     ))
     AMREX_IF_ON_HOST((
         amrex::ignore_unused(offset);
diff --git a/Src/Base/AMReX_IndexType.H b/Src/Base/AMReX_IndexType.H
index 0fd613d2a9..ea310b9b43 100644
--- a/Src/Base/AMReX_IndexType.H
+++ b/Src/Base/AMReX_IndexType.H
@@ -3,145 +3,364 @@
 #define BL_INDEXTYPE_H
 #include <AMReX_Config.H>
 
-#include <AMReX_ccse-mpi.H>
 #include <AMReX_IntVect.H>
 #include <AMReX_SPACE.H>
+#include <AMReX_Tuple.H>
 
 #include <iosfwd>
 
 namespace amrex {
 
+/**
+* \brief Type for defining CellIndex so that all IndexTypeND with different dimensions
+* have the same CellIndex type.
+*/
+struct CellIndexEnum {
+    //! The cell index type: one of CELL or NODE.
+    enum CellIndex { CELL = 0, NODE = 1 };
+};
+
 /**
 * \brief Cell-Based or Node-Based Indices
 *
-* The class IndexType defines an index as being cell based or node (edge)
-* based in each of the AMREX_SPACEDIM directions.  This class defines an
+* The class IndexTypeND defines an index as being cell based or node (edge)
+* based in each of the dim directions. This class defines an
 * enumerated type CellIndex to be either CELL or NODE; i.e. each of the
-* AMREX_SPACEDIM dimensions must be either CELL or NODE.
+* dim dimensions must be either CELL or NODE.
 */
-class IndexType
+template<int dim>
+class IndexTypeND : public CellIndexEnum
 {
-    friend MPI_Datatype ParallelDescriptor::Mpi_typemap<IndexType>::type();
-
 public:
-    //! The cell index type: one of CELL or NODE.
-    enum CellIndex { CELL = 0, NODE = 1 };
+    static_assert(1 <= dim && dim <= 31, "The number of dimensions of IndexTypeND must be positive"
+        " and less than 32");
+
     //! The default constructor
     AMREX_GPU_HOST_DEVICE
-    constexpr IndexType () noexcept  = default;
-    //! Construct an IndexType identical to an IntVect.
+    constexpr IndexTypeND () noexcept = default;
+    //! Construct an IndexTypeND identical to an IntVectND.
     AMREX_GPU_HOST_DEVICE
-    explicit IndexType (const IntVect& iv) noexcept
-        : itype(AMREX_D_TERM((iv[0]?1:0), | ((iv[1]?1:0)<<1), | ((iv[2]?1:0)<<2)))
-        {}
+    explicit IndexTypeND (const IntVectND<dim>& iv) noexcept {
+        for (int i=0; i<dim; ++i) {
+            itype |= (iv[i] ? 1U : 0U) << i;
+        }
+    }
     /**
-    * \brief Construct an IndexType given an explicit CellIndex for
-    * each direction.  AMREX_D_DECL is a macro that sets the constructor
-    * to take AMREX_SPACEDIM arguments.
+    * \brief Construct an IndexTypeND given an explicit CellIndex for
+    * each direction.  The inputs for this constructor are N CellIndex,
+    * where N is equal to the number of dimensions of the IndexTypeND.
     */
+    template <class...Args,
+        std::enable_if_t<
+            (sizeof...(Args)+1 == dim) &&
+            IsConvertible_v<CellIndex, Args...>,
+        int> = 0>
     AMREX_GPU_HOST_DEVICE
-    constexpr IndexType (AMREX_D_DECL(CellIndex i, CellIndex j, CellIndex k)) noexcept
-        : itype(AMREX_D_TERM(i, | (j<<1), | (k<<2)))
-        {}
+    constexpr IndexTypeND (CellIndex i, Args...js) noexcept {
+        CellIndex locarr[dim] = {i, static_cast<CellIndex>(js)...};
+        for (int s=0; s<dim; ++s) {
+            itype |= ((locarr[s] == CellIndex::NODE) ? 1U : 0U) << s;
+        }
+    }
     // dtor, copy-ctor, copy-op=, move-ctor, and move-op= are compiler generated.
 
-    //! Set IndexType to be NODE based in direction dir.
+    //! Set IndexTypeND to be NODE based in direction dir.
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void set (int dir) noexcept { itype |= mask(dir); }
-    //! Set IndexType to be CELL based in direction dir.
+    //! Set IndexTypeND to be CELL based in direction dir.
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void unset (int dir) noexcept { itype &= ~mask(dir); }
-    //! True if IndexType is NODE based in direction dir.
+    //! True if IndexTypeND is NODE based in direction dir.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool test (int dir) const noexcept { return (itype & mask(dir)) != 0; }
     //! Set NODE based in all directions.
-    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    void setall () noexcept { itype = (1 << AMREX_SPACEDIM) - 1; }
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
+    void setall () noexcept { itype = (1U << dim) - 1; }
     //! Set CELL based in all directions.
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void clear () noexcept { itype = 0; }
-    //! True if this IndexType is NODE based in any direction.
+    //! True if this IndexTypeND is NODE based in any direction.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool any () const noexcept { return itype != 0; }
-    //! True if IndexType is valid.
+    //! True if IndexTypeND is valid.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    bool ok () const noexcept { return itype < (1 << AMREX_SPACEDIM); }
+    bool ok () const noexcept { return itype < (1U << dim); }
     //! Change from CELL to NODE or NODE to CELL in direction dir.
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void flip (int i) noexcept { itype ^= mask(i); }
-    //! True if IndexTypes are identical.
+    //! True if IndexTypeNDs are identical.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    bool operator== (const IndexType& t) const noexcept { return t.itype == itype; }
-    //! True if IndexTypes are not identical.
+    bool operator== (const IndexTypeND& t) const noexcept { return t.itype == itype; }
+    //! True if IndexTypeNDs are not identical.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    bool operator!= (const IndexType& t) const noexcept { return t.itype != itype; }
+    bool operator!= (const IndexTypeND& t) const noexcept { return t.itype != itype; }
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    bool operator< (const IndexType& t) const noexcept { return itype < t.itype; }
-    //! True if the IndexType is CELL based in all directions.
+    bool operator< (const IndexTypeND& t) const noexcept { return itype < t.itype; }
+    //! True if the IndexTypeND is CELL based in all directions.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool cellCentered () const noexcept { return itype == 0; }
-    //! True if the IndexType is CELL based in dir-direction.
+    //! True if the IndexTypeND is CELL based in dir-direction.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool cellCentered (int dir) const noexcept { return (itype & mask(dir)) == 0; }
-    //! True if the IndexType is NODE based in all directions.
+    //! True if the IndexTypeND is NODE based in all directions.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    bool nodeCentered () const noexcept { return itype == (1<<AMREX_SPACEDIM)-1; }
-    //! True if the IndexType is NODE based in dir-direction.
+    bool nodeCentered () const noexcept { return itype == (1U<<dim)-1; }
+    //! True if the IndexTypeND is NODE based in dir-direction.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool nodeCentered (int dir) const noexcept { return (itype & mask(dir)) != 0; }
-    //! Set IndexType to CellIndex type t in direction dir.
+    //! Set IndexTypeND to CellIndex type t in direction dir.
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void setType (int dir, CellIndex t) noexcept { t == CELL ? unset(dir) : set(dir); }
     //! Returns  the CellIndex in direction dir.
-    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    CellIndex ixType (int dir) const noexcept { return (CellIndex) ((itype & (1<<dir)) >> dir); }
-    //! Return an integer representing the IndexType in direction dir.
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
+    CellIndex ixType (int dir) const noexcept { return (CellIndex) ((itype & (1U<<dir)) >> dir); }
+    //! Return an integer representing the IndexTypeND in direction dir.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     int operator[] (int dir) const noexcept { return test(dir); }
-    //! Fill an IntVect of size AMREX_SPACEDIM with IndexTypes.
+    //! Returns the i'th CellIndex of the IndexTypeND. Used by structured bindings.
+    template<std::size_t i>
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
+    CellIndex get () const noexcept { static_assert(0<=i && i<dim); return ixType(i); }
+    //! Fill an IntVectND of size dim with IndexTypeNDs.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    IntVect ixType () const noexcept { return IntVect(AMREX_D_DECL(itype&1, (itype>>1)&1, (itype>>2)&1)); }
-    //! Fill an IntVect of size AMREX_SPACEDIM with IndexTypes.
+    IntVectND<dim> ixType () const noexcept {
+        IntVectND<dim> retval(0);
+        for (int i=0; i<dim; ++i) {
+            retval[i] = test(i);
+        }
+        return retval;
+    }
+    //! Fill an IntVectND of size dim with IndexTypeNDs.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    IntVect toIntVect () const noexcept { return IntVect(AMREX_D_DECL(itype&1, (itype>>1)&1, (itype>>2)&1)); }
+    IntVectND<dim> toIntVect () const noexcept {
+        IntVectND<dim> retval(0);
+        for (int i=0; i<dim; ++i) {
+            retval[i] = test(i);
+        }
+        return retval;
+    }
     /**
-    * \brief This static member function returns an IndexType object of value
-    * IndexType::CELL.  It is provided as a convenience to our users
+    * \brief This static member function returns an IndexTypeND object of value
+    * IndexTypeND::CELL.  It is provided as a convenience to our users
     * when defining a Box all of whose faces should be of type
-    * IndexType::CELL.
+    * IndexTypeND::CELL.
     */
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static constexpr IndexType TheCellType () noexcept {
-        return IndexType(AMREX_D_DECL(IndexType::CELL,
-                                      IndexType::CELL,
-                                      IndexType::CELL));
+    static constexpr IndexTypeND<dim> TheCellType () noexcept {
+        return IndexTypeND<dim>{};
     }
     /**
-    * \brief This static member function returns an IndexType object of value
-    * IndexType::NODE.  It is provided as a convenience to our users
+    * \brief This static member function returns an IndexTypeND object of value
+    * IndexTypeND::NODE.  It is provided as a convenience to our users
     * when defining a Box all of whose faces should be of type
-    * IndexType::NODE.
+    * IndexTypeND::NODE.
+    */
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr IndexTypeND<dim> TheNodeType () noexcept {
+        IndexTypeND<dim> retval{};
+        retval.setall();
+        return retval;
+    }
+
+    //! Return the size of this IndexTypeND.
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr std::size_t size () noexcept {
+        return static_cast<std::size_t>(dim);
+    }
+
+    //! Return the size of this IndexTypeND.
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr int isize () noexcept {
+        return dim;
+    }
+
+    using value_type = CellIndex;
+
+    /**
+    * \brief Returns a new IndexTypeND of size new_dim and
+    * assigns the first new_dim values of this IndexTypeND to it.
+    */
+    template<int new_dim>
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    IndexTypeND<new_dim> shrink () const noexcept {
+        static_assert(new_dim <= dim);
+        IndexTypeND<new_dim> retval{};
+        retval.getBits() = itype & ((1U << new_dim) - 1);
+        return retval;
+    }
+
+    /**
+    * \brief Returns a new IndexTypeND of size new_dim and
+    * assigns all values of this IndexTypeND to it and fill_extra to the remaining elements.
+    */
+    template<int new_dim>
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    IndexTypeND<new_dim> expand (CellIndex fill_extra=CellIndex::CELL) const noexcept {
+        static_assert(new_dim >= dim);
+        IndexTypeND<new_dim> retval{};
+        retval.getBits() = itype;
+        if (fill_extra == CellIndex::NODE) {
+            retval.getBits() |= (1U << new_dim) - (1U << dim);
+        }
+        return retval;
+    }
+
+    /**
+    * \brief Returns a new IndexTypeND of size new_dim
+    * by either shrinking or expanding this IndexTypeND
     */
+    template<int new_dim>
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static constexpr IndexType TheNodeType () noexcept {
-        return IndexType(AMREX_D_DECL(IndexType::NODE,
-                                      IndexType::NODE,
-                                      IndexType::NODE));
+    IndexTypeND<new_dim> resize (CellIndex fill_extra=CellIndex::CELL) const noexcept {
+        if constexpr (new_dim > dim) {
+            return expand<new_dim>(fill_extra);
+        } else {
+            return shrink<new_dim>();
+        }
     }
 
+    //! Return the bit field representing the underlying data
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
+    unsigned int& getBits () noexcept { return itype; }
+
+    //! Return the bit field representing the underlying data
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
+    const unsigned int& getBits () const noexcept { return itype; }
+
 private:
     //! Returns 1<<k.
     [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static int mask (int k) noexcept { return 1<<k; }
-    //! An integer holding the CellIndex in bits 0 - AMREX_SPACEDIM-1.
+    static constexpr unsigned int mask (int k) noexcept { return 1U<<k; }
+    //! An integer holding the CellIndex in bits 0 - dim-1.
     unsigned int itype{0};
 };
 
-//! Write an IndexType to an ostream in ASCII.
-std::ostream& operator<< (std::ostream& os, const IndexType& itype);
-//! Read an IndexType from an istream.
-std::istream& operator>> (std::istream& is, IndexType& itype);
+using IndexType = IndexTypeND<AMREX_SPACEDIM>;
+
+// Template deduction guide for IndexTypeND
+template<int dim>
+AMREX_GPU_HOST_DEVICE // __device__ for HIP
+IndexTypeND(const IntVectND<dim>&) -> IndexTypeND<dim>;
+
+// Template deduction guide for IndexTypeND
+template <class...Args,
+    std::enable_if_t<
+        IsConvertible_v<IndexType::CellIndex, Args...>,
+    int> = 0>
+AMREX_GPU_HOST_DEVICE // __device__ for HIP
+IndexTypeND(IndexType::CellIndex, Args...) -> IndexTypeND<sizeof...(Args)+1>;
+
+namespace detail {
+    std::ostream& index_type_write (std::ostream& os, const unsigned int& iv, int dim);
+    std::istream& index_type_read (std::istream& is, unsigned int& iv, int dim);
+
+    template<class T, std::size_t...Ns>
+    [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
+    T IndexTypeSplit_imp (T& retval, std::index_sequence<Ns...>, unsigned int src) noexcept {
+        int dim_shift = 0;
+        (
+            (
+                amrex::get<Ns>(retval).getBits() =
+                    (src >> dim_shift) & ((1U << amrex::get<Ns>(retval).isize()) - 1),
+                dim_shift += amrex::get<Ns>(retval).isize()
+            ), ...
+        );
+        return retval;
+    }
+}
+
+//! Write an IndexTypeND to an ostream in ASCII.
+template<int dim>
+std::ostream& operator<< (std::ostream& os, const IndexTypeND<dim>& it) {
+    return detail::index_type_write(os, it.getBits(), dim);
+}
+//! Read an IndexTypeND from an istream.
+template<int dim>
+std::istream& operator>> (std::istream& is, IndexTypeND<dim>& it) {
+    return detail::index_type_read(is, it.getBits(), dim);
+}
+
+/**
+* \brief Returns a IndexTypeND obtained by concatenating the input IndexTypeNDs.
+* The dimension of the return value equals the sum of the dimensions of the inputted IndexTypeNDs.
+*/
+template<int d, int...dims>
+[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+constexpr IndexTypeND<detail::get_sum<d, dims...>()>
+IndexTypeCat (const IndexTypeND<d>& v, const IndexTypeND<dims>&...vects) noexcept {
+    IndexTypeND<detail::get_sum<d, dims...>()> retval{};
+    retval.getBits() |= v.getBits();
+    int dim_shift = v.isize();
+    (
+        (
+            retval.getBits() |= (vects.getBits() << dim_shift),
+            dim_shift += vects.isize()
+        ), ...
+    );
+    return retval;
+}
+
+/**
+* \brief Returns a tuple of IndexTypeND obtained by splitting the input IndexTypeND
+* according to the dimensions specified by the template arguments.
+*/
+template<int d, int...dims>
+[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+constexpr GpuTuple<IndexTypeND<d>, IndexTypeND<dims>...>
+IndexTypeSplit (const IndexTypeND<detail::get_sum<d, dims...>()>& v) noexcept {
+    GpuTuple<IndexTypeND<d>, IndexTypeND<dims>...> retval{};
+    return detail::IndexTypeSplit_imp(retval,
+                                      std::make_index_sequence<1 + sizeof...(dims)>(),
+                                      v.getBits());
+}
+
+/**
+* \brief Returns a new IndexTypeND of size new_dim and
+* assigns the first new_dim values of v to it.
+*/
+template<int new_dim, int old_dim>
+[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+constexpr IndexTypeND<new_dim>
+IndexTypeShrink (const IndexTypeND<old_dim>& v) noexcept {
+    return v.template shrink<new_dim>();
+}
+
+/**
+* \brief Returns a new IndexTypeND of size new_dim and
+* assigns all values of iv to it and fill_extra to the remaining elements.
+*/
+template<int new_dim, int old_dim>
+[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+constexpr IndexTypeND<new_dim>
+IndexTypeExpand (const IndexTypeND<old_dim>& v,
+                 IndexType::CellIndex fill_extra=IndexType::CellIndex::CELL) noexcept {
+    return v.template expand<new_dim>(fill_extra);
+}
 
+/**
+* \brief Returns a new IndexTypeND of size new_dim
+* by either shrinking or expanding iv.
+*/
+template<int new_dim, int old_dim>
+[[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+constexpr IndexTypeND<new_dim>
+IndexTypeResize (const IndexTypeND<old_dim>& v,
+                 IndexType::CellIndex fill_extra=IndexType::CellIndex::CELL) noexcept {
+    return v.template resize<new_dim>(fill_extra);
 }
 
+} // namespace amrex
+
+// Spcialize std::tuple_size for IndexTypeND. Used by structured bindings.
+template<int dim>
+struct std::tuple_size<amrex::IndexTypeND<dim>> {
+    static constexpr std::size_t value = dim;
+};
+
+// Spcialize std::tuple_element for IndexTypeND. Used by structured bindings.
+template<std::size_t s, int dim>
+struct std::tuple_element<s, amrex::IndexTypeND<dim>> {
+    using type = typename amrex::IndexTypeND<dim>::CellIndex;
+};
+
 #endif /*BL_INDEXTYPE_H*/
diff --git a/Src/Base/AMReX_IndexType.cpp b/Src/Base/AMReX_IndexType.cpp
index e96052a4d1..e2cd7097ec 100644
--- a/Src/Base/AMReX_IndexType.cpp
+++ b/Src/Base/AMReX_IndexType.cpp
@@ -4,16 +4,16 @@
 #include <iostream>
 #include <iomanip>
 
-namespace amrex {
+namespace amrex::detail {
 
 std::ostream&
-operator<< (std::ostream&    os,
-            const IndexType& it)
+index_type_write (std::ostream& os, const unsigned int& iv, int dim)
 {
-    os << '('
-       << AMREX_D_TERM( (it.test(0)?'N':'C'),
-                  << ',' << (it.test(1)?'N':'C'),
-                  << ',' << (it.test(2)?'N':'C')) << ')' << std::flush;
+    os << '(' << (((iv & 1U) != 0) ? 'N' : 'C');
+    for (int i=1; i<dim; ++i) {
+        os << ',' << (((iv & (1U<<i)) != 0) ? 'N' : 'C');
+    }
+    os << ')' << std::flush;
 
     if (os.fail()) {
         amrex::Error("operator<<(ostream&,IndexType&) failed");
@@ -28,19 +28,18 @@ operator<< (std::ostream&    os,
 #define BL_IGNORE_MAX 100000
 
 std::istream&
-operator>> (std::istream& is,
-            IndexType&    it)
+index_type_read (std::istream& is, unsigned int& iv, int dim)
 {
-    char AMREX_D_DECL(t0,t1,t2);
-
-    AMREX_D_EXPR( is.ignore(BL_IGNORE_MAX, '(') >> t0,
-            is.ignore(BL_IGNORE_MAX, ',') >> t1,
-            is.ignore(BL_IGNORE_MAX, ',') >> t2);
+    char t = '0';
+    is.ignore(BL_IGNORE_MAX, '(') >> t;
+    BL_ASSERT(t == 'C' || t == 'N');
+    t == 'N' ? (iv |= 1U) : (iv &= ~1U);
+    for (int i=1; i<dim; ++i) {
+        is.ignore(BL_IGNORE_MAX, ',') >> t;
+        BL_ASSERT(t == 'C' || t == 'N');
+        t == 'N' ? (iv |= (1U << i)) : (iv &= ~(1U << i));
+    }
     is.ignore(BL_IGNORE_MAX, ')');
-    AMREX_D_TERM(
-        BL_ASSERT(t0 == 'C' || t0 == 'N'); t0=='N'?it.set(0):it.unset(0); ,
-        BL_ASSERT(t1 == 'C' || t1 == 'N'); t1=='N'?it.set(1):it.unset(1); ,
-        BL_ASSERT(t2 == 'C' || t2 == 'N'); t2=='N'?it.set(2):it.unset(2));
 
     if (is.fail()) {
         amrex::Error("operator>>(ostream&,IndexType&) failed");
diff --git a/Src/Base/AMReX_IntVect.H b/Src/Base/AMReX_IntVect.H
index 2cf48053cd..c133a1f57c 100644
--- a/Src/Base/AMReX_IntVect.H
+++ b/Src/Base/AMReX_IntVect.H
@@ -1108,7 +1108,7 @@ namespace detail {
     template<class T, std::size_t...Ns>
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr
     T IntVectSplit_imp (T& retval, std::index_sequence<Ns...>, const int * src) noexcept {
-        (IntVectSplit_imp2(std::get<Ns>(retval), src), ...);
+        (IntVectSplit_imp2(amrex::get<Ns>(retval), src), ...);
         return retval;
     }
 
diff --git a/Src/Base/AMReX_MultiFabUtil.H b/Src/Base/AMReX_MultiFabUtil.H
index d4b3a3271c..228070a13c 100644
--- a/Src/Base/AMReX_MultiFabUtil.H
+++ b/Src/Base/AMReX_MultiFabUtil.H
@@ -398,6 +398,20 @@ namespace amrex
      * \param stddev standard deviation of normal distribution
      */
     void FillRandomNormal (MultiFab& mf, int scomp, int ncomp, Real mean, Real stddev);
+
+    /**
+     * \brief Convexify AMR data
+     *
+     * This function "convexifies" the AMR data by removing cells that are
+     * covered by fine levels from coarse level MultiFabs. This could be
+     * useful for visualization. The returned MultiFabs on coarse levels
+     * have different BoxArrays from the original BoxArrays. For the finest
+     * level, the data is simply copied to the returned object. The returned
+     * MultiFabs have no ghost cells. For nodal data, the nodes on the
+     * coarse/fine interface exist on both levels.
+     */
+    [[nodiscard]] Vector<MultiFab> convexify (Vector<MultiFab const*> const& mf,
+                                Vector<IntVect> const& refinement_ratio);
 }
 
 namespace amrex {
diff --git a/Src/Base/AMReX_MultiFabUtil.cpp b/Src/Base/AMReX_MultiFabUtil.cpp
index a2a1f14ed8..86a1e29054 100644
--- a/Src/Base/AMReX_MultiFabUtil.cpp
+++ b/Src/Base/AMReX_MultiFabUtil.cpp
@@ -1216,4 +1216,57 @@ namespace amrex
             FillRandomNormal(p, npts, mean, stddev);
         }
     }
+
+    Vector<MultiFab> convexify (Vector<MultiFab const*> const& mf,
+                                Vector<IntVect> const& refinement_ratio)
+    {
+        if (mf.empty()) { return Vector<MultiFab>{}; }
+
+        const auto nlevels = int(mf.size());
+        Vector <MultiFab> rmf(nlevels);
+
+        const int ncomp = mf[nlevels-1]->nComp();
+        rmf[nlevels-1].define(mf[nlevels-1]->boxArray(),
+                              mf[nlevels-1]->DistributionMap(), ncomp, 0);
+        MultiFab::Copy(rmf[nlevels-1], *mf[nlevels-1], 0, 0, ncomp, 0);
+
+        for (int ilev = nlevels-2; ilev >= 0; --ilev) {
+            BoxArray fba = mf[ilev+1]->boxArray();
+            BoxArray cba = mf[ilev  ]->boxArray();
+            AMREX_ASSERT(fba.ixType() == cba.ixType());
+            AMREX_ASSERT(mf[ilev]->nComp() == ncomp);
+
+            fba.convert(IntVect(0)).coarsen(refinement_ratio[ilev]);
+            cba.convert(IntVect(0));
+            auto const& cdm = mf[ilev]->DistributionMap();
+
+            BoxList blnew, bltmp;
+            Vector<int> procmap;
+            Vector<int> localmap;
+            for (int ibox = 0; ibox < int(cba.size()); ++ibox) {
+                fba.complementIn(bltmp, cba[ibox]);
+                blnew.join(bltmp);
+                procmap.resize(procmap.size()+bltmp.size(), cdm[ibox]);
+                if (ParallelDescriptor::MyProc() == cdm[ibox]) {
+                    localmap.resize(localmap.size()+bltmp.size(), ibox);
+                }
+            }
+
+            if (blnew.isNotEmpty()) {
+                BoxArray banew(std::move(blnew));
+                banew.convert(mf[ilev]->ixType());
+                DistributionMapping dmnew(std::move(procmap));
+                rmf[ilev].define(banew, dmnew, ncomp, 0);
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+                for (MFIter mfi(rmf[ilev], TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+                   rmf[ilev][mfi].template copy<RunOn::Device>
+                       ((*mf[ilev])[localmap[mfi.LocalIndex()]], mfi.tilebox());
+                }
+            }
+        }
+
+        return rmf;
+    }
 }
diff --git a/Src/Base/AMReX_PODVector.H b/Src/Base/AMReX_PODVector.H
index 0e7458fe5b..86791b09c0 100644
--- a/Src/Base/AMReX_PODVector.H
+++ b/Src/Base/AMReX_PODVector.H
@@ -366,7 +366,7 @@ namespace amrex
             else
             {
                 // if the allocators are not the same we give up and copy
-                *this = a_vector;
+                *this = a_vector; // must copy instead of move
             }
 
             return *this;
diff --git a/Src/Base/AMReX_PlotFileDataImpl.cpp b/Src/Base/AMReX_PlotFileDataImpl.cpp
index b85c17ad93..28cc38e27f 100644
--- a/Src/Base/AMReX_PlotFileDataImpl.cpp
+++ b/Src/Base/AMReX_PlotFileDataImpl.cpp
@@ -31,7 +31,7 @@ PlotFileDataImpl::PlotFileDataImpl (std::string const& plotfile_name)
     for (int i = 0; i < m_ncomp; ++i) {
         std::string tmp;
         std::getline(is, tmp);
-        m_var_names[i] = amrex::trim(tmp);
+        m_var_names[i] = amrex::trim(std::move(tmp));
     }
 
     is >> m_spacedim >> m_time >> m_finest_level;
diff --git a/Src/Base/AMReX_Random.cpp b/Src/Base/AMReX_Random.cpp
index a0b72155c5..891a69e140 100644
--- a/Src/Base/AMReX_Random.cpp
+++ b/Src/Base/AMReX_Random.cpp
@@ -187,15 +187,14 @@ UniqueRandomSubset (Vector<int> &uSet, int setSize, int poolSize,
     Abort("**** Error in UniqueRandomSubset:  setSize > poolSize.");
   }
   std::set<int> copySet;
-  Vector<int> uSetTemp;
+  uSet.clear();
   while(static_cast<int>(copySet.size()) < setSize) {
     int r = static_cast<int>(Random_int(poolSize));
     if(copySet.find(r) == copySet.end()) {
       copySet.insert(r);
-      uSetTemp.push_back(r);
+      uSet.push_back(r);
     }
   }
-  uSet = uSetTemp;
   if(printSet) {
     for(int i(0); i < uSet.size(); ++i) {
         AllPrint() << "uSet[" << i << "]  = " << uSet[i] << '\n';
diff --git a/Src/EB/AMReX_EB2_GeometryShop.H b/Src/EB/AMReX_EB2_GeometryShop.H
index 1be28d536c..1278931005 100644
--- a/Src/EB/AMReX_EB2_GeometryShop.H
+++ b/Src/EB/AMReX_EB2_GeometryShop.H
@@ -241,7 +241,7 @@ public:
         {
             const auto& problo = geom.ProbLoArray();
             const auto& dx = geom.CellSizeArray();
-            auto f = m_f;
+            auto const& f = m_f;
             ReduceOps<ReduceOpSum,ReduceOpSum,ReduceOpSum> reduce_op;
             ReduceData<int,int,int> reduce_data(reduce_op);
             using ReduceTuple = typename decltype(reduce_data)::Type;
@@ -297,13 +297,13 @@ public:
     void fillFab (BaseFab<Real>& levelset, const Geometry& geom, RunOn run_on,
                   Box const& bounding_box) const noexcept
     {
-        const auto problo = geom.ProbLoArray();
-        const auto dx = geom.CellSizeArray();
+        const auto& problo = geom.ProbLoArray();
+        const auto& dx = geom.CellSizeArray();
         const Box& bx = levelset.box();
         const auto& a = levelset.array();
-        const auto blo = amrex::lbound(bounding_box);
-        const auto bhi = amrex::ubound(bounding_box);
-        auto f = m_f;
+        const auto& blo = amrex::lbound(bounding_box);
+        const auto& bhi = amrex::ubound(bounding_box);
+        const auto& f = m_f;
         AMREX_HOST_DEVICE_FOR_3D_FLAG(run_on, bx, i, j, k,
         {
             a(i,j,k) = f(AMREX_D_DECL(problo[0]+amrex::Clamp(i,blo.x,bhi.x)*dx[0],
@@ -335,11 +335,11 @@ public:
     void fillFab_Cpu (BaseFab<Real>& levelset, const Geometry& geom,
                       Box const& bounding_box) const noexcept
     {
-        const auto problo = geom.ProbLoArray();
-        const auto dx = geom.CellSizeArray();
+        const auto& problo = geom.ProbLoArray();
+        const auto& dx = geom.CellSizeArray();
         const Box& bx = levelset.box();
-        const auto blo = amrex::lbound(bounding_box);
-        const auto bhi = amrex::ubound(bounding_box);
+        const auto& blo = amrex::lbound(bounding_box);
+        const auto& bhi = amrex::ubound(bounding_box);
 
         const auto& a = levelset.array();
         amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
@@ -359,9 +359,9 @@ public:
     {
         auto const& dx = geom.CellSizeArray();
         auto const& problo = geom.ProbLoArray();
-        const auto blo = amrex::lbound(bounding_box);
-        const auto bhi = amrex::ubound(bounding_box);
-        auto f = m_f;
+        const auto& blo = amrex::lbound(bounding_box);
+        const auto& bhi = amrex::ubound(bounding_box);
+        auto const& f = m_f;
         for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
             Array4<Real> const& inter = inter_arr[idim];
             Array4<Type_t const> const& type = type_arr[idim];
@@ -430,8 +430,8 @@ public:
     {
         auto const& dx = geom.CellSizeArray();
         auto const& problo = geom.ProbLoArray();
-        const auto blo = amrex::lbound(bounding_box);
-        const auto bhi = amrex::ubound(bounding_box);
+        const auto& blo = amrex::lbound(bounding_box);
+        const auto& bhi = amrex::ubound(bounding_box);
         const Box bx{inter};
         amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
         {
diff --git a/Src/EB/AMReX_EBDataCollection.cpp b/Src/EB/AMReX_EBDataCollection.cpp
index 141f061462..f9b3d66c79 100644
--- a/Src/EB/AMReX_EBDataCollection.cpp
+++ b/Src/EB/AMReX_EBDataCollection.cpp
@@ -208,8 +208,8 @@ void EBDataCollection::extendDataOutsideDomain (IntVect const& level_ng)
                     if (apbx.bigEnd(idim) == nbx.bigEnd(idim)) {
                         apbx.growHi(idim,-1);
                     }
-                    auto lev_apidim_domain = lev_ap_domain[idim];
-                    Dim3 off = IntVect::TheDimensionVector(idim).dim3();
+                    auto const& lev_apidim_domain = lev_ap_domain[idim];
+                    Dim3 const& off = IntVect::TheDimensionVector(idim).dim3();
                     amrex::ParallelFor(apbx,
                     [=] AMREX_GPU_DEVICE (int i, int j, int k)
                     {
diff --git a/Src/EB/AMReX_distFcnElement.H b/Src/EB/AMReX_distFcnElement.H
index 1d7c2ef059..d0c7b86713 100644
--- a/Src/EB/AMReX_distFcnElement.H
+++ b/Src/EB/AMReX_distFcnElement.H
@@ -18,11 +18,11 @@ class distFcnElement2d { // NOLINT(cppcoreguidelines-special-member-functions)
 
   [[nodiscard]] virtual amrex::Real cpdist(amrex::RealVect pt, amrex::RealVect& cp) const = 0;
   [[nodiscard]] virtual amrex::Real cpside(amrex::RealVect pt, amrex::RealVect& cp) const = 0;
-  static int solve_thomas(const std::vector<amrex::Real> &a,
-                   std::vector<amrex::Real> b,
-                   const std::vector<amrex::Real> &c,
-                   std::vector<amrex::Real> d,
-                   std::vector<amrex::Real> &X);
+  static int solve_thomas (const std::vector<amrex::Real> &a,
+                           std::vector<amrex::Real> &b,
+                           const std::vector<amrex::Real> &c,
+                           std::vector<amrex::Real> &d,
+                           std::vector<amrex::Real> &X);
 };
 
 
diff --git a/Src/EB/AMReX_distFcnElement.cpp b/Src/EB/AMReX_distFcnElement.cpp
index b6040d30f5..c72bb6af8f 100644
--- a/Src/EB/AMReX_distFcnElement.cpp
+++ b/Src/EB/AMReX_distFcnElement.cpp
@@ -6,11 +6,11 @@
 
 namespace amrex {
 
-int distFcnElement2d::solve_thomas(const std::vector<amrex::Real> &a,
-                                   std::vector<amrex::Real> b,
-                                   const std::vector<amrex::Real> &c,
-                                   std::vector<amrex::Real> d,
-                                   std::vector<amrex::Real> &x)
+int distFcnElement2d::solve_thomas (const std::vector<amrex::Real> &a,
+                                    std::vector<amrex::Real> &b,
+                                    const std::vector<amrex::Real> &c,
+                                    std::vector<amrex::Real> &d,
+                                    std::vector<amrex::Real> &x)
 {
   int n = static_cast<int>(d.size());
   x.resize(n);
diff --git a/Src/Extern/HYPRE/AMReX_HypreABecLap.cpp b/Src/Extern/HYPRE/AMReX_HypreABecLap.cpp
index 4607fcc1fc..2d4c313c08 100644
--- a/Src/Extern/HYPRE/AMReX_HypreABecLap.cpp
+++ b/Src/Extern/HYPRE/AMReX_HypreABecLap.cpp
@@ -115,7 +115,7 @@ HypreABecLap::getSolution (MultiFab& a_soln)
         auto reghi = Hypre::hiV(reg);
         HYPRE_StructVectorGetBoxValues(x, reglo.data(), reghi.data(), (*soln)[mfi].dataPtr());
     }
-    Gpu::synchronize();
+    Gpu::hypreSynchronize();
 
     if (a_soln.nGrowVect() != 0) {
         MultiFab::Copy(a_soln, tmp, 0, 0, 1, 0);
@@ -235,7 +235,7 @@ HypreABecLap::prepareSolver ()
         HYPRE_StructMatrixSetBoxValues(A, reglo.data(), reghi.data(),
                                        regular_stencil_size, stencil_indices.data(),
                                        mat);
-        Gpu::synchronize();
+        Gpu::hypreSynchronize();
     }
     HYPRE_StructMatrixAssemble(A);
 
@@ -299,7 +299,7 @@ HypreABecLap::loadVectors (MultiFab& soln, const MultiFab& rhs)
         HYPRE_StructVectorSetBoxValues(x, reglo.data(), reghi.data(), soln[mfi].dataPtr());
         HYPRE_StructVectorSetBoxValues(b, reglo.data(), reghi.data(), rhs_diag[mfi].dataPtr());
     }
-    Gpu::synchronize();
+    Gpu::hypreSynchronize();
 }
 
 }
diff --git a/Src/Extern/HYPRE/AMReX_HypreABecLap2.cpp b/Src/Extern/HYPRE/AMReX_HypreABecLap2.cpp
index 1b7ba8e895..a44d2d517e 100644
--- a/Src/Extern/HYPRE/AMReX_HypreABecLap2.cpp
+++ b/Src/Extern/HYPRE/AMReX_HypreABecLap2.cpp
@@ -137,7 +137,7 @@ HypreABecLap2::getSolution (MultiFab& a_soln)
         HYPRE_SStructVectorGetBoxValues(x, part, reglo.data(), reghi.data(),
                                         0, (*soln)[mfi].dataPtr());
     }
-    Gpu::synchronize();
+    Gpu::hypreSynchronize();
 
     if (a_soln.nGrowVect() != 0) {
         MultiFab::Copy(a_soln, tmp, 0, 0, 1, 0);
@@ -262,7 +262,7 @@ HypreABecLap2::prepareSolver ()
         HYPRE_SStructMatrixSetBoxValues(A, part, reglo.data(), reghi.data(),
                                         0, regular_stencil_size, stencil_indices.data(),
                                         mat);
-        Gpu::synchronize();
+        Gpu::hypreSynchronize();
     }
     HYPRE_SStructMatrixAssemble(A);
 
@@ -335,7 +335,7 @@ HypreABecLap2::loadVectors (MultiFab& soln, const MultiFab& rhs)
         HYPRE_SStructVectorSetBoxValues(b, part, reglo.data(), reghi.data(),
                                         0, rhs_diag[mfi].dataPtr());
     }
-    Gpu::synchronize();
+    Gpu::hypreSynchronize();
 }
 
 }
diff --git a/Src/Extern/HYPRE/AMReX_HypreABecLap3.cpp b/Src/Extern/HYPRE/AMReX_HypreABecLap3.cpp
index d04b85c314..216cf368f5 100644
--- a/Src/Extern/HYPRE/AMReX_HypreABecLap3.cpp
+++ b/Src/Extern/HYPRE/AMReX_HypreABecLap3.cpp
@@ -77,7 +77,7 @@ HypreABecLap3::getSolution (MultiFab& a_soln)
             (*l_soln)[mfi].setVal<RunOn::Device>(0.0);
         }
     }
-    Gpu::synchronize();
+    Gpu::hypreSynchronize();
 
     if (use_tmp_mf) {
         MultiFab::Copy(a_soln, tmp, 0, 0, 1, 0);
@@ -499,9 +499,9 @@ HypreABecLap3::prepareSolver ()
                 });
             }
 
-            Gpu::synchronize();
+            Gpu::streamSynchronize();
             HYPRE_IJMatrixSetValues(A,nrows,ncols,rows,cols,mat);
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
         }
     }
     HYPRE_IJMatrixAssemble(A);
@@ -681,7 +681,7 @@ HypreABecLap3::loadVectors (MultiFab& soln, const MultiFab& rhs)
             HYPRE_IJVectorSetValues(b, nrows, cell_id_vec[mfi].dataPtr(), rhs_diag[mfi].dataPtr());
         }
     }
-    Gpu::synchronize();
+    Gpu::hypreSynchronize();
 }
 
 }  // namespace amrex
diff --git a/Src/Extern/HYPRE/AMReX_HypreNodeLap.cpp b/Src/Extern/HYPRE/AMReX_HypreNodeLap.cpp
index bb4bc557ca..33e45ca066 100644
--- a/Src/Extern/HYPRE/AMReX_HypreNodeLap.cpp
+++ b/Src/Extern/HYPRE/AMReX_HypreNodeLap.cpp
@@ -123,9 +123,9 @@ HypreNodeLap::HypreNodeLap (const BoxArray& grids_, const DistributionMapping& d
                 adjust_singular_matrix(ncols, cols, rows, mat);
             }
 
-            Gpu::synchronize();
+            Gpu::streamSynchronize();
             HYPRE_IJMatrixSetValues(A, nrows, ncols, rows, cols, mat);
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
         }
     }
     HYPRE_IJMatrixAssemble(A);
@@ -324,9 +324,9 @@ HypreNodeLap::loadVectors (MultiFab& soln, const MultiFab& rhs)
                 });
             }
 
-            Gpu::synchronize();
+            Gpu::streamSynchronize();
             HYPRE_IJVectorSetValues(b, nrows, rows_vec.data(), bvec.data());
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
         }
     }
 }
@@ -347,7 +347,7 @@ HypreNodeLap::getSolution (MultiFab& soln)
             xvec.resize(nrows);
             Real* xp = xvec.data();
             HYPRE_IJVectorGetValues(x, nrows, rows_vec.data(), xp);
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
 
             const Box& bx = mfi.validbox();
             const auto& xfab = tmpsoln.array(mfi);
@@ -359,7 +359,7 @@ HypreNodeLap::getSolution (MultiFab& soln)
                 }
             });
 
-            Gpu::synchronize();
+            Gpu::streamSynchronize();
         }
     }
 
diff --git a/Src/Extern/HYPRE/AMReX_HypreSolver.H b/Src/Extern/HYPRE/AMReX_HypreSolver.H
index 8f34b72f2d..fb70c90bbb 100644
--- a/Src/Extern/HYPRE/AMReX_HypreSolver.H
+++ b/Src/Extern/HYPRE/AMReX_HypreSolver.H
@@ -569,7 +569,7 @@ HypreSolver<MSS>::fill_matrix (Filler const& filler)
             Gpu::streamSynchronize();
             HYPRE_IJMatrixSetValues(m_A, nrows, ncols_vec.data(), rows,
                                     cols_vec.data(), mat_vec.data());
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
         }
     }
     HYPRE_IJMatrixAssemble(m_A);
@@ -653,7 +653,7 @@ HypreSolver<MSS>::load_vectors (Vector<MF      *> const& a_soln,
             Gpu::streamSynchronize();
             HYPRE_IJVectorSetValues(m_x, nrows, rows, xp);
             HYPRE_IJVectorSetValues(m_b, nrows, rows, bp);
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
         }
     }
 }
@@ -683,7 +683,7 @@ HypreSolver<MSS>::get_solution (Vector<MF*> const& a_soln)
             HYPRE_Int const* rows = m_global_id_vec[mfi].data();
 
             HYPRE_IJVectorGetValues(m_x, nrows, rows, xp);
-            Gpu::synchronize();
+            Gpu::hypreSynchronize();
 
             HYPRE_Int offset = 0;
             for (int ivar = 0; ivar < m_nvars; ++ivar) {
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H
index db8a10f208..b613a4f3a8 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H
@@ -396,6 +396,7 @@ MLCGSolverT<MF>::solve_cg (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
         if ( !initial_vec_zeroed ) {
             LocalAdd(sol, sorig, 0, 0, ncomp, nghost);
         }
+        if (ret == 8) { ret = 9; }
     }
     else
     {
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H
index 02f7adeb11..970cf48fc2 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H
@@ -562,8 +562,8 @@ MLCellABecLapT<MF>::addInhomogNeumannFlux (
             if (! domain.contains(ccb)) {
                 for (int icomp = 0; icomp < ncomp; ++icomp) {
                     auto const& phi = sol.const_array(mfi,icomp);
-                    auto const bv = bndry.bndryValues(ori).multiFab().const_array(mfi,icomp);
-                    auto const bc = bcoef[idim] ? bcoef[idim]->const_array(mfi,icomp)
+                    auto const& bv = bndry.bndryValues(ori).multiFab().const_array(mfi,icomp);
+                    auto const& bc = bcoef[idim] ? bcoef[idim]->const_array(mfi,icomp)
                         : Array4<RT const>{};
                     auto const& f = grad[idim]->array(mfi,icomp);
                     if (ori.isLow()) {
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H
index 77f9abe409..c02bde6d04 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H
@@ -1558,29 +1558,32 @@ MLMGT<MF>::actualBottomSolve ()
             } else {
                 cg_type = MLCGSolverT<MF>::Type::BiCGStab;
             }
+
             int ret = bottomSolveWithCG(x, *bottom_b, cg_type);
-            // If the MLMG solve failed then set the correction to zero
-            if (ret != 0 && ret != 9) {
+
+            if (ret != 0 && (bottom_solver == BottomSolver::cgbicg ||
+                             bottom_solver == BottomSolver::bicgcg))
+            {
+                if (bottom_solver == BottomSolver::cgbicg) {
+                    cg_type = MLCGSolverT<MF>::Type::BiCGStab; // switch to bicg
+                } else {
+                    cg_type = MLCGSolverT<MF>::Type::CG; // switch to cg
+                }
                 setVal(cor[amrlev][mglev], RT(0.0));
-                if (bottom_solver == BottomSolver::cgbicg ||
-                    bottom_solver == BottomSolver::bicgcg) {
-                    if (bottom_solver == BottomSolver::cgbicg) {
-                        cg_type = MLCGSolverT<MF>::Type::BiCGStab; // switch to bicg
+                ret = bottomSolveWithCG(x, *bottom_b, cg_type);
+                if (ret == 0) { // switch permanently
+                    if (cg_type == MLCGSolverT<MF>::Type::CG) {
+                        bottom_solver = BottomSolver::cg;
                     } else {
-                        cg_type = MLCGSolverT<MF>::Type::CG; // switch to cg
-                    }
-                    ret = bottomSolveWithCG(x, *bottom_b, cg_type);
-                    if (ret != 0) {
-                        setVal(cor[amrlev][mglev], RT(0.0));
-                    } else { // switch permanently
-                        if (cg_type == MLCGSolverT<MF>::Type::CG) {
-                            bottom_solver = BottomSolver::cg;
-                        } else {
-                            bottom_solver = BottomSolver::bicgstab;
-                        }
+                        bottom_solver = BottomSolver::bicgstab;
                     }
                 }
             }
+
+            // If the bottom solve failed then set the correction to zero
+            if (ret != 0 && ret != 9) {
+                setVal(cor[amrlev][mglev], RT(0.0));
+            }
             const int n = (ret==0) ? nub : nuf;
             for (int i = 0; i < n; ++i) {
                 linop.smooth(amrlev, mglev, x, b);
diff --git a/Src/Particle/AMReX_ParticleContainerI.H b/Src/Particle/AMReX_ParticleContainerI.H
index cd2d030e55..de76a104b5 100644
--- a/Src/Particle/AMReX_ParticleContainerI.H
+++ b/Src/Particle/AMReX_ParticleContainerI.H
@@ -71,9 +71,9 @@ IntVect
 ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>::Index (const P& p, int lev) const
 {
     const Geometry& geom = Geom(lev);
-    const auto domain = geom.Domain();
-    const auto plo = geom.ProbLoArray();
-    const auto dxi = geom.InvCellSizeArray();
+    const auto& domain = geom.Domain();
+    const auto& plo = geom.ProbLoArray();
+    const auto& dxi = geom.InvCellSizeArray();
 
     return Assignor{}(p, plo, dxi, domain);
 }
diff --git a/Src/Particle/AMReX_ParticleUtil.H b/Src/Particle/AMReX_ParticleUtil.H
index e40d43c9b3..6b31db7e99 100644
--- a/Src/Particle/AMReX_ParticleUtil.H
+++ b/Src/Particle/AMReX_ParticleUtil.H
@@ -58,9 +58,9 @@ numParticlesOutOfRange (Iterator const& pti, IntVect nGrow)
     const auto& ptd = tile.getConstParticleTileData();
     const auto& geom = pti.Geom(pti.GetLevel());
 
-    const auto domain = geom.Domain();
-    const auto plo = geom.ProbLoArray();
-    const auto dxi = geom.InvCellSizeArray();
+    const auto& domain = geom.Domain();
+    const auto& plo = geom.ProbLoArray();
+    const auto& dxi = geom.InvCellSizeArray();
 
     Box box = pti.tilebox();
     box.grow(nGrow);