Porting fortran files to cpp in Src/Extern/HYPRE (AMReX-Codes#1049)

Src/Extern/HYPRE/AMReX_Habec_K.H

@@ -0,0 +1,122 @@
+#ifndef AMREX_Habec_K_H_
+#define AMREX_Habec_K_H_
+
+#if (AMREX_SPACEDIM == 2)
+#include <AMReX_Habec_2D_K.H>
+#elif (AMREX_SPACEDIM == 3)
+#include <AMReX_Habec_3D_K.H>
+#endif
+
+namespace amrex {
+
+AMREX_FORCE_INLINE
+void amrex_hpacoef (Box const& box,
+                    BaseFab<GpuArray<Real,2*AMREX_SPACEDIM + 1>>& mat,
+                    const FArrayBox& a,
+                    Real& sa)
+{
+    Array4<GpuArray<Real,2*AMREX_SPACEDIM + 1>> const& mat_arr = mat.array();
+    Array4<Real const> const& a_arr = a.const_array();
+    if (sa == 0.0){
+        AMREX_HOST_DEVICE_PARALLEL_FOR_3D (box, i, j, k,
+        {
+            mat_arr(i,j,k)[0] = 0.0;
+        });
+    }else{
+        AMREX_HOST_DEVICE_PARALLEL_FOR_3D (box, i, j, k,
+        {
+            mat_arr(i,j,k)[0] = sa * a_arr(i,j,k);
+        });
+    }
+}
+
+AMREX_FORCE_INLINE
+void amrex_hpdiag (Box const& box,
+                   BaseFab<GpuArray<Real,2*AMREX_SPACEDIM + 1>>& mat,
+                   FArrayBox& diag)
+{
+    Array4<GpuArray<Real,2*AMREX_SPACEDIM + 1>> const& mat_arr = mat.array();
+    Array4<Real> const& diag_arr = diag.array();
+
+    AMREX_HOST_DEVICE_PARALLEL_FOR_3D (box, i, j, k,
+    {
+        diag_arr(i,j,k) = 1.0/mat_arr(i,j,k)[0];
+        for (int ii=0; ii<2*AMREX_SPACEDIM + 1; ii++){
+            mat_arr(i,j,k)[ii] = mat_arr(i,j,k)[ii]*diag_arr(i,j,k);
+        }
+    });
+}
+
+#ifdef AMREX_USE_EB
+
+template<typename T> AMREX_FORCE_INLINE
+void amrex_hpeb_fill_cellid (Box const& box,
+                             T& nrows,
+                             BaseFab<T>& cell_id,
+                             const EBCellFlagFab& flag)
+{
+    Array4<T> const& cell_id_arr = cell_id.array();
+    Array4<const EBCellFlag> const& flag_arr = flag.array();
+
+    nrows = 0;
+    Gpu::DeviceScalar<int> nrows_gpu(nrows);
+    int* nrowsg = nrows_gpu.dataPtr();
+
+    AMREX_HOST_DEVICE_PARALLEL_FOR_3D (box, i, j, k,
+    {
+        if (!flag_arr(i,j,k).isCovered()){
+            cell_id_arr(i,j,k) = *nrowsg;
+            Gpu::Atomic::Add(nrowsg, 1);
+        }
+    });
+    nrows = nrows_gpu.dataValue();
+}
+
+AMREX_FORCE_INLINE
+void amrex_hpeb_copy_from_vec (Box const& box,
+                               FArrayBox& a,
+                               Real* v,
+                               const EBCellFlagFab& flag)
+{
+    Array4<Real> const& a_arr = a.array();
+    Array4<const EBCellFlag> const& flag_arr = flag.array();
+
+    int nrows = 0;
+    Gpu::DeviceScalar<int> nrows_gpu(nrows);
+    int* nrowsg = nrows_gpu.dataPtr();
+
+    AMREX_HOST_DEVICE_PARALLEL_FOR_3D (box, i, j, k,
+    {
+        if (!flag_arr(i,j,k).isCovered()){
+            a_arr(i,j,k) = v[*nrowsg];
+            Gpu::Atomic::Add(nrowsg, 1);
+        }
+    });
+}
+
+AMREX_FORCE_INLINE
+void amrex_hpeb_copy_to_vec (Box const& box,
+                             FArrayBox& a,
+                             Real* v,
+                             const EBCellFlagFab& flag)
+{
+    Array4<Real> const& a_arr = a.array();
+    Array4<const EBCellFlag> const& flag_arr = flag.array();
+
+    int nrows = 0;
+    Gpu::DeviceScalar<int> nrows_gpu(nrows);
+    int* nrowsg = nrows_gpu.dataPtr();
+
+    AMREX_HOST_DEVICE_PARALLEL_FOR_3D (box, i, j, k,
+    {
+        if (!flag_arr(i,j,k).isCovered()){
+            v[*nrowsg] = a_arr(i,j,k);
+            Gpu::Atomic::Add(nrowsg, 1);
+        }
+    });
+}
+
+#endif
+
+}
+#endif

Src/Extern/HYPRE/AMReX_HypreABecLap.cpp

@@ -1,9 +1,10 @@
 
 #include <AMReX_HypreABecLap.H>
-#include <AMReX_HypreABec_F.H>
 #include <string>
 #include <algorithm>
 
+#include <AMReX_Habec_K.H>
+
 #include <_hypre_struct_mv.h>
 
 namespace amrex {
@@ -195,24 +196,16 @@ HypreABecLap::prepareSolver ()
     std::iota(stencil_indices.begin(), stencil_indices.end(), 0);
     const Real* dx = geom.CellSize();
     const int bho = (m_maxorder > 2) ? 1 : 0;
-    FArrayBox rfab;
+    BaseFab<GpuArray<Real,regular_stencil_size>> rfab;
     for (MFIter mfi(acoefs); mfi.isValid(); ++mfi)
     {  
         const Box &reg = mfi.validbox();
 
-        rfab.resize(reg,regular_stencil_size);
-        Real* mat = rfab.dataPtr();
-
-        amrex_hpacoef(BL_TO_FORTRAN_BOX(reg),
-                      mat,
-                      BL_TO_FORTRAN_ANYD(acoefs[mfi]),
-                      &scalar_a);
+        rfab.resize(reg);
+        amrex_hpacoef(reg, rfab, acoefs[mfi], scalar_a); 
 
         for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
-            amrex_hpbcoef(BL_TO_FORTRAN_BOX(reg),
-                          mat,
-                          BL_TO_FORTRAN_ANYD(bcoefs[idim][mfi]),
-                          &scalar_b, dx, &idim);
+            amrex_hpbcoef(reg, rfab, bcoefs[idim][mfi], scalar_b, dx, idim);
         }
 
         const Vector< Vector<BoundCond> > & bcs_i = m_bndry->bndryConds(mfi);
@@ -227,17 +220,12 @@ HypreABecLap::prepareSolver ()
             const int bctype = bcs_i[cdir][0];
             const Real  &bcl = bcl_i[cdir];
             const Mask  &msk = m_bndry->bndryMasks(ori)[mfi];
-
-            amrex_hpmat(BL_TO_FORTRAN_BOX(reg),
-                        mat,
-                        BL_TO_FORTRAN_ANYD(bcoefs[idim][mfi]),
-                        BL_TO_FORTRAN_ANYD(msk),
-                        &scalar_b, dx, &cdir, &bctype, &bcl, &bho);
+
+            amrex_hpmat(reg, rfab, bcoefs[idim][mfi], msk, scalar_b, dx, cdir, bctype, bcl, bho);   
         }
-
-        amrex_hpdiag(BL_TO_FORTRAN_BOX(reg),
-                     mat,
-                     BL_TO_FORTRAN_ANYD(diaginv[mfi]));
+
+        amrex_hpdiag(reg, rfab, diaginv[mfi]);
+        Real* mat = (Real*) rfab.dataPtr();
 
         auto reglo = Hypre::loV(reg);
         auto reghi = Hypre::hiV(reg);

Src/Extern/HYPRE/AMReX_HypreABecLap2.cpp

@@ -1,5 +1,6 @@
 #include <AMReX_HypreABecLap2.H>
-#include <AMReX_HypreABec_F.H>
+
+#include <AMReX_Habec_K.H>
 
 #include <cmath>
 #include <numeric>
@@ -219,24 +220,16 @@ HypreABecLap2::prepareSolver ()
     const HYPRE_Int part = 0;
     const Real* dx = geom.CellSize();
     const int bho = (m_maxorder > 2) ? 1 : 0;
-    FArrayBox rfab;
+    BaseFab<GpuArray<Real, regular_stencil_size>> rfab;
     for (MFIter mfi(acoefs); mfi.isValid(); ++mfi)
     {
         const Box &reg = mfi.validbox();
 
-        rfab.resize(reg,regular_stencil_size);
-        Real* mat = rfab.dataPtr();
-
-        amrex_hpacoef(BL_TO_FORTRAN_BOX(reg),
-                      mat,
-                      BL_TO_FORTRAN_ANYD(acoefs[mfi]),
-                      &scalar_a);
+        rfab.resize(reg);
+        amrex_hpacoef(reg, rfab, acoefs[mfi], scalar_a);
 
         for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
-            amrex_hpbcoef(BL_TO_FORTRAN_BOX(reg),
-                          mat,
-                          BL_TO_FORTRAN_ANYD(bcoefs[idim][mfi]),
-                          &scalar_b, dx, &idim);
+            amrex_hpbcoef(reg, rfab, bcoefs[idim][mfi], scalar_b, dx, idim);
         }
 
         const Vector< Vector<BoundCond> > & bcs_i = m_bndry->bndryConds(mfi);
@@ -251,16 +244,11 @@ HypreABecLap2::prepareSolver ()
             const Real &bcl  = bcl_i[cdir];
             const Mask &msk  = m_bndry->bndryMasks(ori)[mfi];
 
-            amrex_hpmat(BL_TO_FORTRAN_BOX(reg),
-                        mat,
-                        BL_TO_FORTRAN_ANYD(bcoefs[idim][mfi]),
-                        BL_TO_FORTRAN_ANYD(msk),
-                        &scalar_b, dx, &cdir, &bctype, &bcl, &bho);
+            amrex_hpmat(reg, rfab, bcoefs[idim][mfi], msk, scalar_b, dx, cdir, bctype, bcl, bho);
         }
 
-        amrex_hpdiag(BL_TO_FORTRAN_BOX(reg),
-                     mat,
-                     BL_TO_FORTRAN_ANYD(diaginv[mfi]));
+        amrex_hpdiag(reg, rfab, diaginv[mfi]); 
+        Real* mat = (Real*) rfab.dataPtr();
 
         // initialize matrix
         auto reglo = Hypre::loV(reg);

Src/Extern/HYPRE/AMReX_HypreABecLap3.H

@@ -38,7 +38,7 @@ private :
     HYPRE_Solver solver = NULL;
 
     LayoutData<HYPRE_Int> ncells_grid;
-    LayoutData<Vector<HYPRE_Int> > cell_id_vec;
+    LayoutData<Gpu::ManagedDeviceVector<HYPRE_Int> > cell_id_vec;
     FabArray<BaseFab<HYPRE_Int> > cell_id;
 
     MultiFab const* m_eb_b_coeffs = nullptr;