Improve AMReX diagnostics and kernels

swm_amrex/CMakeLists.txt

@@ -31,7 +31,7 @@ find_package(AMReX REQUIRED ${AMREX_COMPONENTS})
 add_subdirectory(swm_AMReX)
 #add_subdirectory(swm_AMReX_Fkernels)
 
-# These version are only for offloading onto the GPU with CUDA backend
-if(SWM_DEVICE STREQUAL "gpu" AND SWM_CUDA)
+# These versions are only for offloading onto the GPU with CUDA backend
+if(SWM_DEVICE STREQUAL "gpu" AND SWM_CUDA AND SWM_FORTRAN)
   add_subdirectory(swm_AMReX_Fsubroutine)
-endif()
+endif()

swm_amrex/common/main.cpp

@@ -60,6 +60,13 @@ int main (int argc, char* argv[])
 
     amrex::MultiFab p;
     DefineNodalMultiFab(psi, p);
+
+    amrex::Print() << "\nGrid Summary:\n";
+    PrintGridSummary(psi, "psi (cell-centered)");
+    PrintGridSummary(p, "p (nodal)        ");
+    PrintGridSummary(u, "u (y-face)       ");
+    PrintGridSummary(v, "v (x-face)       ");
+    amrex::Print() << std::endl;
 
     // **********************************
     // Initialize Data
@@ -71,6 +78,8 @@ int main (int argc, char* argv[])
 
     InitializeVariables(geom, psi, p, u, v);
 
+    WriteDiagonalElements("diagonal_initial.dat", p, u, v, nx);
+
     // **********************************
     // Write initial plot file
     // **********************************
@@ -178,6 +187,8 @@ int main (int argc, char* argv[])
 
     BL_PROFILE_VAR_STOP(total);
 
+    WriteDiagonalElements("diagonal_final.dat", p, u, v, nx);
+
     }
 
     amrex::Finalize();

swm_amrex/common/swm_mini_app_utils.cpp

@@ -1,7 +1,12 @@
 #include <string>
 #include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <limits>
+#include <vector>
 
 #include <AMReX.H>
+#include <AMReX_Gpu.H>
 #include <AMReX_ParmParse.H>
 #include <AMReX_PlotFileUtil.H>
 #include <AMReX_MultiFab.H>
@@ -221,7 +226,7 @@ void InitializeVariables(const amrex::Geometry & geom,
             const amrex::Real x_transformed = LinearMapCoordinates(x_node, x_min, x_max, 0.0, 2*pi);
             const amrex::Real y_transformed = LinearMapCoordinates(y_node, y_min, y_max, 0.0, 2*pi);
 
-            p_array(i,j,k) = pcf * (std::cos(2*x_transformed) + std::cos(2*y_transformed)) + 5000;
+            p_array(i,j,k) = pcf * (std::cos(2*x_transformed) + std::cos(2*y_transformed)) + 50000.;
         });
     }
 
@@ -349,3 +354,102 @@ void UpdateVariables(const amrex::Geometry& geom,
 
     return;
 }
+
+void PrintGridSummary(const amrex::MultiFab& mf, const std::string& name)
+{
+    const amrex::BoxArray& ba = mf.boxArray();
+
+    int numgrid = static_cast<int>(ba.size());
+    amrex::Long ncells = ba.numPts();
+    amrex::Box minbox = ba.minimalBox();
+    double ntot = static_cast<double>(minbox.numPts());
+    amrex::Real frac = amrex::Real(100.0 * double(ncells) / ntot);
+
+    // AMReX-style output format
+    amrex::Print() << "  " << name
+                   << "   "
+                   << numgrid
+                   << " grids  "
+                   << ncells
+                   << " cells  "
+                   << frac
+                   << " % of domain  "
+                   << mf.nGrow()
+                   << " ghost"
+                   << std::endl;
+
+    if (numgrid > 1) {
+        amrex::Long vmin = std::numeric_limits<amrex::Long>::max();
+        amrex::Long vmax = -1;
+        for (int i = 0; i < numgrid; ++i) {
+            amrex::Long v = ba[i].numPts();
+            vmin = std::min(vmin, v);
+            vmax = std::max(vmax, v);
+        }
+        amrex::Print() << "            "
+                       << "smallest grid: " << vmin
+                       << "  largest grid: " << vmax
+                       << std::endl;
+    }
+}
+
+void WriteDiagonalElements(const std::string& filename,
+                           const amrex::MultiFab& p,
+                           const amrex::MultiFab& u,
+                           const amrex::MultiFab& v,
+                           int n_diag)
+{
+    // Get the valid box to determine array bounds
+    const amrex::Box& p_box = p.boxArray().minimalBox();
+    int nx = p_box.length(0);
+    int ny = p_box.length(1);
+    int mnmin = std::min(std::min(nx, ny), n_diag);
+
+    // Create host-readable copies of the MultiFabs.
+    amrex::MultiFab p_host(p.boxArray(), p.DistributionMap(), p.nComp(), 0,
+                           amrex::MFInfo().SetArena(amrex::The_Pinned_Arena()));
+    amrex::MultiFab u_host(u.boxArray(), u.DistributionMap(), u.nComp(), 0,
+                           amrex::MFInfo().SetArena(amrex::The_Pinned_Arena()));
+    amrex::MultiFab v_host(v.boxArray(), v.DistributionMap(), v.nComp(), 0,
+                           amrex::MFInfo().SetArena(amrex::The_Pinned_Arena()));
+
+    amrex::MultiFab::Copy(p_host, p, 0, 0, p.nComp(), 0);
+    amrex::MultiFab::Copy(u_host, u, 0, 0, u.nComp(), 0);
+    amrex::MultiFab::Copy(v_host, v, 0, 0, v.nComp(), 0);
+    amrex::Gpu::streamSynchronize();
+
+    // Collect diagonal values
+    std::vector<amrex::Real> p_diag(mnmin), u_diag(mnmin), v_diag(mnmin);
+    for (int i = 0; i < mnmin; i++) {
+        for (amrex::MFIter mfi(p_host); mfi.isValid(); ++mfi) {
+            if (mfi.validbox().contains(amrex::IntVect(i, i))) {
+                p_diag[i] = p_host.const_array(mfi)(i, i, 0);
+            }
+        }
+        for (amrex::MFIter mfi(u_host); mfi.isValid(); ++mfi) {
+            if (mfi.validbox().contains(amrex::IntVect(i, i))) {
+                u_diag[i] = u_host.const_array(mfi)(i, i, 0);
+            }
+        }
+        for (amrex::MFIter mfi(v_host); mfi.isValid(); ++mfi) {
+            if (mfi.validbox().contains(amrex::IntVect(i, i))) {
+                v_diag[i] = v_host.const_array(mfi)(i, i, 0);
+            }
+        }
+    }
+
+    // Write to file
+    std::ofstream ofs(filename);
+    if (ofs) {
+        ofs << std::setprecision(15);
+        ofs << "p:";
+        for (int i = 0; i < mnmin; i++) ofs << " " << p_diag[i];
+        ofs << "\nu:";
+        for (int i = 0; i < mnmin; i++) ofs << " " << u_diag[i];
+        ofs << "\nv:";
+        for (int i = 0; i < mnmin; i++) ofs << " " << v_diag[i];
+        ofs << "\n";
+        ofs.close();
+        amrex::Print() << " diagonal elements written to " << filename << std::endl;
+    }
+}

swm_amrex/common/swm_mini_app_utils.h

@@ -1,6 +1,8 @@
 #ifndef SWM_MINI_APP_UTILS_H_
 #define SWM_MINI_APP_UTILS_H_
 
+#include <string>
+
 #include <AMReX.H>
 #include <AMReX_MultiFab.H>
 
@@ -77,4 +79,12 @@ void UpdateOldVariables(const double alpha, const int time_step, const amrex::Ge
                         const amrex::MultiFab& p_new, const amrex::MultiFab& u_new, const amrex::MultiFab& v_new, 
                         amrex::MultiFab& p_old, amrex::MultiFab& u_old, amrex::MultiFab& v_old);
 
+void WriteDiagonalElements(const std::string& filename,
+                           const amrex::MultiFab& p,
+                           const amrex::MultiFab& u,
+                           const amrex::MultiFab& v,
+                           int n_diag);
+
+void PrintGridSummary(const amrex::MultiFab& mf, const std::string& name);
+
 #endif // SWM_MINI_APP_UTILS_H_

swm_amrex/doc/figure_labels.tex

@@ -97,7 +97,7 @@ \section{Initial Condition}
 
 
 The initial pressure evaluated at all the nodal points using:
-$$ P_0 =  \frac{\pi^2 A^2}{L_x^2} \left[ \cos(2 \xi_x) + \cos(2 \xi_y) \right] + 5000 $$ 
+$$ P_0 =  \frac{\pi^2 A^2}{L_x^2} \left[ \cos(2 \xi_x) + \cos(2 \xi_y) \right] + 50000 $$
 
 
 
@@ -458,4 +458,4 @@ \subsection{Equations for interpolating from one location to another}
 	\end{minipage}
 \end{center}
 
-\end{document}
+\end{document}

swm_amrex/swm_AMReX/swm_mini_app_kernels.h

@@ -14,10 +14,20 @@ void UpdateIntermediateVariablesKernel( const int i, const int j, const int k,
                                         const amrex::Array4<amrex::Real>& h,
                                         const amrex::Array4<amrex::Real>& z)
 {
-    cu(i,j,k) = 0.5*(p(i,j,k) + p(i+1,j,k))*u(i,j,k);
-    cv(i,j,k) = 0.5*(p(i,j,k) + p(i,j+1,k))*v(i,j,k);
-    z(i,j,k) = (fsdx*(v(i+1,j,k)-v(i,j,k)) - fsdy*(u(i,j+1,k)-u(i,j,k)))/(p(i,j,k)+p(i+1,j,k)+p(i,j+1,k)+p(i+1,j+1,k));
-    h(i,j,k) = p(i,j,k) + 0.25*(u(i-1,j,k)*u(i-1,j,k) + u(i,j,k)*u(i,j,k) + v(i,j-1,k)*v(i,j-1,k) + v(i,j,k)*v(i,j,k));
+    amrex::Real p_ijk = p(i,j,k);
+    amrex::Real p_i1jk = p(i+1,j,k);
+    amrex::Real p_ij1k = p(i,j+1,k);
+    amrex::Real p_i1j1k = p(i+1,j+1,k);
+    amrex::Real u_ijk = u(i,j,k);
+    amrex::Real u_im1jk = u(i-1,j,k);
+    amrex::Real u_ij1k = u(i,j+1,k);
+    amrex::Real v_ijk = v(i,j,k);
+    amrex::Real v_i1jk = v(i+1,j,k);
+    amrex::Real v_ijm1k = v(i,j-1,k);
+    cu(i,j,k) = 0.5*(p_ijk + p_i1jk)*u_ijk;
+    cv(i,j,k) = 0.5*(p_ijk + p_ij1k)*v_ijk;
+    z(i,j,k) = (fsdx*(v_i1jk - v_ijk) - fsdy*(u_ij1k - u_ijk))/(p_ijk + p_i1jk + p_ij1k + p_i1j1k);
+    h(i,j,k) = p_ijk + 0.25*(u_im1jk*u_im1jk + u_ijk*u_ijk + v_ijm1k*v_ijm1k + v_ijk*v_ijk);
 }
 
 
@@ -35,9 +45,23 @@ void UpdateNewVariablesKernel( const int i, const int j, const int k,
                                const amrex::Array4<amrex::Real>& u_new,
                                const amrex::Array4<amrex::Real>& v_new)
 {
-    u_new(i,j,k) = u_old(i,j,k) + tdts8 * (z(i,j-1,k)+z(i,j,k)) * (cv(i,j-1,k) + cv(i,j,k) + cv(i+1,j-1,k) + cv(i+1,j,k)) - tdtsdx * (h(i+1,j,k) - h(i,j,k));
-    v_new(i,j,k) = v_old(i,j,k) - tdts8 * (z(i-1,j,k)+z(i,j,k)) * (cu(i-1,j,k) + cu(i-1,j+1,k) + cu(i,j,k) + cu(i,j+1,k)) - tdtsdy * (h(i,j+1,k) - h(i,j,k));
-    p_new(i,j,k) = p_old(i,j,k) - tdtsdx * (cu(i,j,k) - cu(i-1,j,k)) - tdtsdy * (cv(i,j,k) - cv(i,j-1,k));
+    amrex::Real z_ijk = z(i,j,k);
+    amrex::Real z_ijm1k = z(i,j-1,k);
+    amrex::Real z_im1jk = z(i-1,j,k);
+    amrex::Real cv_ijk = cv(i,j,k);
+    amrex::Real cv_ijm1k = cv(i,j-1,k);
+    amrex::Real cv_i1jk = cv(i+1,j,k);
+    amrex::Real cv_i1jm1k = cv(i+1,j-1,k);
+    amrex::Real cu_ijk = cu(i,j,k);
+    amrex::Real cu_im1jk = cu(i-1,j,k);
+    amrex::Real cu_ij1k = cu(i,j+1,k);
+    amrex::Real cu_im1j1k = cu(i-1,j+1,k);
+    amrex::Real h_ijk = h(i,j,k);
+    amrex::Real h_i1jk = h(i+1,j,k);
+    amrex::Real h_ij1k = h(i,j+1,k);
+    u_new(i,j,k) = u_old(i,j,k) + tdts8 * (z_ijm1k + z_ijk) * (cv_ijm1k + cv_ijk + cv_i1jm1k + cv_i1jk) - tdtsdx * (h_i1jk - h_ijk);
+    v_new(i,j,k) = v_old(i,j,k) - tdts8 * (z_im1jk + z_ijk) * (cu_im1jk + cu_im1j1k + cu_ijk + cu_ij1k) - tdtsdy * (h_ij1k - h_ijk);
+    p_new(i,j,k) = p_old(i,j,k) - tdtsdx * (cu_ijk - cu_im1jk) - tdtsdy * (cv_ijk - cv_ijm1k);
 }
 
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE