diff --git a/swm_amrex/CMakeLists.txt b/swm_amrex/CMakeLists.txt index 6ba97ba..357939c 100644 --- a/swm_amrex/CMakeLists.txt +++ b/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() \ No newline at end of file +endif() diff --git a/swm_amrex/common/main.cpp b/swm_amrex/common/main.cpp index f78d5d3..a2de7c9 100644 --- a/swm_amrex/common/main.cpp +++ b/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(); diff --git a/swm_amrex/common/swm_mini_app_utils.cpp b/swm_amrex/common/swm_mini_app_utils.cpp index 50a1d11..bdaeb57 100644 --- a/swm_amrex/common/swm_mini_app_utils.cpp +++ b/swm_amrex/common/swm_mini_app_utils.cpp @@ -1,7 +1,12 @@ #include #include +#include +#include +#include +#include #include +#include #include #include #include @@ -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(ba.size()); + amrex::Long ncells = ba.numPts(); + amrex::Box minbox = ba.minimalBox(); + double ntot = static_cast(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::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 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; + } +} diff --git a/swm_amrex/common/swm_mini_app_utils.h b/swm_amrex/common/swm_mini_app_utils.h index be3e791..f1e6529 100644 --- a/swm_amrex/common/swm_mini_app_utils.h +++ b/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 + #include #include @@ -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_ diff --git a/swm_amrex/doc/figure_labels.tex b/swm_amrex/doc/figure_labels.tex index c51358f..f9b834b 100644 --- a/swm_amrex/doc/figure_labels.tex +++ b/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} \ No newline at end of file +\end{document} diff --git a/swm_amrex/swm_AMReX/swm_mini_app_kernels.h b/swm_amrex/swm_AMReX/swm_mini_app_kernels.h index bbc1088..f30da60 100644 --- a/swm_amrex/swm_AMReX/swm_mini_app_kernels.h +++ b/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& h, const amrex::Array4& 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& u_new, const amrex::Array4& 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