Merge branch 'development' into coriolis

hklion · Jan 5, 2024 · 1f98b86 · 1f98b86
2 parents 9533e6f + a3c5017
commit 1f98b86
Show file tree

Hide file tree

Showing 50 changed files with 3,194 additions and 4,397 deletions.
diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst
@@ -71,6 +71,9 @@ Notes
 
 -  The NeTCDF option is only available if ROMSX has been built with USE_NETCDF enabled.
 
+-  Velocity components are defined on faces within the ROMS-X code, but are averaged onto
+   cell centers when written in amrex/native plotfiles.
+
 -  File prefixes can include directories.
 
 .. _examples-of-usage-8:

diff --git a/Exec/OCCAMS/prob.cpp b/Exec/OCCAMS/prob.cpp
@@ -82,18 +82,18 @@ init_custom_coriolis (const Geometry& geom,
 
 void
 init_custom_prob(
-        const Box& bx,
-        Array4<Real      > const& state,
-        Array4<Real      > const& x_vel,
-        Array4<Real      > const& y_vel,
-        Array4<Real      > const& z_vel,
+        const Box& /*bx*/,
+        Array4<Real      > const& /*state*/,
+        Array4<Real      > const& /*x_vel*/,
+        Array4<Real      > const& /*y_vel*/,
+        Array4<Real      > const& /*z_vel*/,
         Array4<Real const> const& /*z_w*/,
-        Array4<Real const> const& z_r,
+        Array4<Real const> const& /*z_r*/,
         Array4<Real const> const& /*Hz*/,
         Array4<Real const> const& /*h*/,
         Array4<Real const> const& /*Zt_avg1*/,
-        GeometryData const& geomdata,
-        const SolverChoice& m_solverChoice)
+        GeometryData const& /*geomdata*/,
+        const SolverChoice& /*m_solverChoice*/)
 {
     Abort("Shouldn't be in init_custom_prob!");
 }

diff --git a/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp b/Source/BoundaryConditions/BoundaryConditions_netcdf.cpp
@@ -74,52 +74,46 @@ ROMSX::fill_from_bdyfiles (MultiFab& mf_to_fill, const Real time, const int bdy_
         {
             Box mf_box(mf_to_fill[mfi.index()].box());
 
-            // We only proceed if the ghost cells of this FAB are not contained by domain
-            if (!domain.contains(mf_box)) {
-
-                // Compute intersections of the FAB to be filled and the bdry data boxes
-                Box xlo = bdy_data_xlo[n_time][ivar].box() & mf_box;
-                Box xhi = bdy_data_xhi[n_time][ivar].box() & mf_box;
-                Box ylo = bdy_data_ylo[n_time][ivar].box() & mf_box;
-                Box yhi = bdy_data_yhi[n_time][ivar].box() & mf_box;
-
-                const Array4<Real>& dest_arr = mf_to_fill.array(mfi,icomp);
-
-                if (!xlo.isEmpty()) {
-                    ParallelFor(xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                    {
-                        dest_arr(i,j,k,icomp) = oma   * bdatxlo_n  (i,j,k,0)
-                                              + alpha * bdatxlo_np1(i,j,k,0);
-                    });
-                }
-
-                if (!xhi.isEmpty()) {
-                    ParallelFor(xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                    {
-                        dest_arr(i,j,k,icomp) = oma   * bdatxhi_n  (i,j,k,0)
-                                              + alpha * bdatxhi_np1(i,j,k,0);
-                    });
-                }
-
-                if (!ylo.isEmpty()) {
-                    ParallelFor(ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                    {
-                        dest_arr(i,j,k,icomp) = oma   * bdatylo_n  (i,j,k,0)
-                                              + alpha * bdatylo_np1(i,j,k,0);
-                    });
-                }
-
-                if (!yhi.isEmpty()) {
-                    ParallelFor(yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                    {
-                        dest_arr(i,j,k,icomp) = oma    * bdatyhi_n  (i,j,k,0)
-                                               + alpha * bdatyhi_np1(i,j,k,0);
-                    });
-                }
-
-            } // FAB not contained in domain
+            // Compute intersections of the FAB to be filled and the bdry data boxes
+            Box xlo = bdy_data_xlo[n_time][ivar].box() & mf_box;
+            Box xhi = bdy_data_xhi[n_time][ivar].box() & mf_box;
+            Box ylo = bdy_data_ylo[n_time][ivar].box() & mf_box;
+            Box yhi = bdy_data_yhi[n_time][ivar].box() & mf_box;
+
+            const Array4<Real>& dest_arr = mf_to_fill.array(mfi,icomp);
+
+            if (!xlo.isEmpty()) {
+                ParallelFor(xlo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    dest_arr(i,j,k,icomp) = oma   * bdatxlo_n  (i,j,k,0)
+                                          + alpha * bdatxlo_np1(i,j,k,0);
+                });
+            }
+
+            if (!xhi.isEmpty()) {
+                ParallelFor(xhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    dest_arr(i,j,k,icomp) = oma   * bdatxhi_n  (i,j,k,0)
+                                          + alpha * bdatxhi_np1(i,j,k,0);
+                });
+            }
+
+            if (!ylo.isEmpty()) {
+                ParallelFor(ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    dest_arr(i,j,k,icomp) = oma   * bdatylo_n  (i,j,k,0)
+                                          + alpha * bdatylo_np1(i,j,k,0);
+                });
+            }
+
+            if (!yhi.isEmpty()) {
+                ParallelFor(yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                {
+                    dest_arr(i,j,k,icomp) = oma    * bdatyhi_n  (i,j,k,0)
+                                           + alpha * bdatyhi_np1(i,j,k,0);
+                });
+            }
         } // mfi
     } // icomp
-
 }
 #endif
diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp
@@ -10,7 +10,7 @@ using namespace amrex;
 //     so this follows the BCVars enum
 //
 void ROMSXPhysBCFunct::impose_xvel_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
-                                        const GpuArray<Real,AMREX_SPACEDIM> dxInv,
+                                        const GpuArray<Real,AMREX_SPACEDIM> /*dxInv*/,
                                         Real /*time*/, int bccomp)
 {
     BL_PROFILE_VAR("impose_xvel_bcs()",impose_xvel_bcs);

diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp
@@ -10,7 +10,7 @@ using namespace amrex;
 //     so this follows the BCVars enum
 //
 void ROMSXPhysBCFunct::impose_yvel_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
-                                        const GpuArray<Real,AMREX_SPACEDIM> dxInv,
+                                        const GpuArray<Real,AMREX_SPACEDIM> /*dxInv*/,
                                         Real /*time*/, int bccomp)
 {
     BL_PROFILE_VAR("impose_yvel_bcs()",impose_yvel_bcs);

diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp
@@ -11,10 +11,8 @@ using namespace amrex;
 //     so this follows the BCVars enum
 //
 void ROMSXPhysBCFunct::impose_zvel_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
-                                        const Array4<Real const>& velx_arr,
-                                        const Array4<Real const>& vely_arr,
-                                        const GpuArray<Real,AMREX_SPACEDIM> dx,
-                                        const GpuArray<Real,AMREX_SPACEDIM> dxInv,
+                                        const GpuArray<Real,AMREX_SPACEDIM> /*dx*/,
+                                        const GpuArray<Real,AMREX_SPACEDIM> /*dxInv*/,
                                         Real /*time*/, int bccomp)
 {
     const auto& dom_lo = amrex::lbound(domain);

diff --git a/Source/BoundaryConditions/ROMSX_FillPatch.cpp b/Source/BoundaryConditions/ROMSX_FillPatch.cpp
@@ -14,7 +14,11 @@ PhysBCFunctNoOp null_bc;
 //
 void
 ROMSX::FillPatch (int lev, Real time, MultiFab& mf_to_fill, Vector<MultiFab*> const& mfs,
+#ifdef ROMSX_USE_NETCDF
                   const int bdy_var_type)
+#else
+                  const int /*bdy_var_type*/)
+#endif
 {
     BL_PROFILE_VAR("ROMSX::FillPatch()",ROMSX_FillPatch);
     int bccomp;

diff --git a/Source/BoundaryConditions/ROMSX_PhysBCFunct.H b/Source/BoundaryConditions/ROMSX_PhysBCFunct.H
@@ -71,8 +71,6 @@ public:
 
     void impose_zvel_bcs (const amrex::Array4<amrex::Real>& dest_arr,
                           const amrex::Box& bx, const amrex::Box& domain,
-                          const amrex::Array4<amrex::Real const>& velx_arr,
-                          const amrex::Array4<amrex::Real const>& vely_arr,
                           const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> dx,
                           const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> dxInv,
                           amrex::Real time, int bccomp);

diff --git a/Source/BoundaryConditions/ROMSX_PhysBCFunct.cpp b/Source/BoundaryConditions/ROMSX_PhysBCFunct.cpp
@@ -24,7 +24,6 @@ void ROMSXPhysBCFunct::operator() (MultiFab& mf, int icomp, int ncomp, IntVect c
     BL_PROFILE("ROMSXPhysBCFunct::()");
 
     const auto& domain = m_geom.Domain();
-    const auto dx      = m_geom.CellSizeArray();
     const auto dxInv   = m_geom.InvCellSizeArray();
 
     // Create a grown domain box containing valid + periodic cells

diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
@@ -138,31 +138,13 @@ ROMSX::WritePlotFile (int which, Vector<std::string> plot_var_names)
             }
         }
 
-        // Next, check for velocities and if desired, output them
+        // Next, check for velocities and output none or all
         if (containerHasElement(plot_var_names, "x_velocity") ||
             containerHasElement(plot_var_names, "y_velocity") ||
             containerHasElement(plot_var_names, "z_velocity")) {
-            amrex::Print()<<"For now, print faces as if they are at cell centers"<<std::endl;
-            //            average_face_to_cellcenter(mf[lev],mf_comp,
-            //                Array<const MultiFab*,3>{&xvel_new[lev],&yvel_new[lev],&zvel_new[lev]});
-            //
-            // Convert the map-factor-scaled-velocities back to velocities
-            //
-            MultiFab dmf(mf[lev], make_alias, mf_comp, AMREX_SPACEDIM);
-            for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-            {
-                const Box& bx = mfi.growntilebox();
-                const Array4<Real> vel_arr = dmf.array(mfi);
-                const Array4<const Real> velx_arr = xvel_new[lev]->const_array(mfi);
-                const Array4<const Real> vely_arr = yvel_new[lev]->const_array(mfi);
-                const Array4<const Real> velz_arr = zvel_new[lev]->const_array(mfi);
-                ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
-                {
-                    vel_arr(i,j,k,0) = velx_arr(i,j,k);
-                    vel_arr(i,j,k,1) = vely_arr(i,j,k);
-                    vel_arr(i,j,k,2) = velz_arr(i,j,k);
-                });
-            }
+            amrex::Print()<<"We now average the face-based velocity components onto cell centers for plotting "<<std::endl;
+            average_face_to_cellcenter(mf[lev],mf_comp,
+                                       Array<const MultiFab*,3>{xvel_new[lev],yvel_new[lev],zvel_new[lev]});
             mf_comp += AMREX_SPACEDIM;
         }
 
@@ -184,10 +166,10 @@ ROMSX::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
                     loc_arr(i,j,k,0) = (i+0.5) * dx;
                     loc_arr(i,j,k,1) = (j+0.5) * dy;
-                    loc_arr(i,j,k,2) = Real(0.125) * (zp_arr(i,j  ,k  ) + zp_arr(i+1,j  ,k  ) +
-                                                      zp_arr(i,j+1,k  ) + zp_arr(i+1,j+1,k  ) +
-                                                      zp_arr(i,j  ,k+1) + zp_arr(i+1,j  ,k+1) +
-                                                      zp_arr(i,j+1,k+1) + zp_arr(i+1,j+1,k+1) );
+                    loc_arr(i,j,k,2) = 0.125_rt * (zp_arr(i,j  ,k  ) + zp_arr(i+1,j  ,k  ) +
+                                                   zp_arr(i,j+1,k  ) + zp_arr(i+1,j+1,k  ) +
+                                                   zp_arr(i,j  ,k+1) + zp_arr(i+1,j  ,k+1) +
+                                                   zp_arr(i,j+1,k+1) + zp_arr(i+1,j+1,k+1) );
                 });
             } // mfi
             mf_comp += AMREX_SPACEDIM;
@@ -250,8 +232,8 @@ ROMSX::WritePlotFile (int which, Vector<std::string> plot_var_names)
 #ifdef ROMSX_USE_NETCDF
         } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") {
              int lev   = 0;
-             int which = 0;
-             writeNCPlotFile(lev, which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]);
+             int nc_which = 0;
+             writeNCPlotFile(lev, nc_which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]);
              total_plot_file_step_1 += 1;
 #endif
         } else {
@@ -318,8 +300,8 @@ ROMSX::WritePlotFile (int which, Vector<std::string> plot_var_names)
 #ifdef ROMSX_USE_NETCDF
         } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") {
              for (int lev = 0; lev <= finest_level; ++lev) {
-                 for (int which = 0; which < num_boxes_at_level[lev]; which++) {
-                     writeNCPlotFile(lev, which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]);
+                 for (int nc_which = 0; nc_which < num_boxes_at_level[lev]; nc_which++) {
+                     writeNCPlotFile(lev, nc_which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]);
                      total_plot_file_step_1 += 1;
                  }
              }

diff --git a/Source/IO/ReadFromInitNetcdf.cpp b/Source/IO/ReadFromInitNetcdf.cpp
@@ -6,7 +6,7 @@ using namespace amrex;
 
 #ifdef ROMSX_USE_NETCDF
 void
-read_data_from_netcdf (int lev,
+read_data_from_netcdf (int /*lev*/,
                        const Box& domain,
                        const std::string& fname,
                        FArrayBox& NC_temp_fab, FArrayBox& NC_salt_fab,
@@ -33,7 +33,7 @@ read_data_from_netcdf (int lev,
 }
 
 void
-read_bathymetry_from_netcdf (int lev,
+read_bathymetry_from_netcdf (int /*lev*/,
                              const Box& domain,
                              const std::string& fname,
                              FArrayBox& NC_h_fab,

diff --git a/Source/Initialization/ROMSX_init_from_netcdf.cpp b/Source/Initialization/ROMSX_init_from_netcdf.cpp
@@ -13,7 +13,7 @@ using namespace amrex;
 #ifdef ROMSX_USE_NETCDF
 
 void
-read_data_from_netcdf (int lev, const Box& domain, const std::string& fname,
+read_data_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
                        FArrayBox& NC_temp_fab, FArrayBox& NC_salt_fab,
                        FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab,
                        FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab,
@@ -125,8 +125,6 @@ ROMSX::init_bathymetry_from_netcdf (int lev)
     Vector<FArrayBox> NC_pm_fab    ; NC_pm_fab.resize(num_boxes_at_level[lev]);
     Vector<FArrayBox> NC_pn_fab    ; NC_pn_fab.resize(num_boxes_at_level[lev]);
 
-    int nboxes = NC_h_fab.size();
-
     for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
     {
         read_bathymetry_from_netcdf(lev, boxes_at_level[lev][idx], nc_grid_file[lev][idx],
@@ -208,7 +206,7 @@ ROMSX::init_coriolis_from_netcdf (int lev)
  * @param lev Integer specifying the current level
  */
 void
-ROMSX::init_bdry_from_netcdf (int lev)
+ROMSX::init_bdry_from_netcdf ()
 {
     if (nc_bdry_file.empty()) {
         amrex::Error("NetCDF boundary file name must be provided via input");

diff --git a/Source/ROMSX.H b/Source/ROMSX.H
@@ -690,7 +690,7 @@ private:
     int bdy_set_width{0};
 
     void init_data_from_netcdf (int lev);
-    void init_bdry_from_netcdf (int lev);
+    void init_bdry_from_netcdf ();
     void init_bathymetry_from_netcdf (int lev);
     void init_coriolis_from_netcdf (int lev);
 #endif // ROMSX_USE_NETCDF

diff --git a/Source/ROMSX.cpp b/Source/ROMSX.cpp
@@ -429,7 +429,7 @@ ROMSX::init_only(int lev, Real time)
         amrex::Print() << "Initial data loaded from netcdf file \n " << std::endl;
 
         amrex::Print() << "Calling init_bdry_from_netcdf " << std::endl;
-        init_bdry_from_netcdf(lev);
+        init_bdry_from_netcdf();
         amrex::Print() << "Boundary data loaded from netcdf file \n " << std::endl;
 #endif
     } else {

diff --git a/Source/ROMSX_make_new_level.cpp b/Source/ROMSX_make_new_level.cpp
@@ -302,20 +302,20 @@ void ROMSX::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& d
     set_weights(lev);
 
     //consider tracking ru and rv indexes more specifically or more similarly to indx
-    vec_ru[lev]->setVal(0.0);
-    vec_rv[lev]->setVal(0.0);
-
-    vec_DU_avg1[lev]->setVal(0.0);
-    vec_DU_avg2[lev]->setVal(0.0);
-    vec_DV_avg1[lev]->setVal(0.0);
-    vec_DV_avg2[lev]->setVal(0.0);
-    vec_rubar[lev]->setVal(0.0);
-    vec_rvbar[lev]->setVal(0.0);
-    vec_rzeta[lev]->setVal(0.0);
-
-    vec_ubar[lev]->setVal(0.0);
-    vec_vbar[lev]->setVal(0.0);
-    vec_zeta[lev]->setVal(0.0);
+    vec_ru[lev]->setVal(0.0_rt);
+    vec_rv[lev]->setVal(0.0_rt);
+
+    vec_DU_avg1[lev]->setVal(0.0_rt);
+    vec_DU_avg2[lev]->setVal(0.0_rt);
+    vec_DV_avg1[lev]->setVal(0.0_rt);
+    vec_DV_avg2[lev]->setVal(0.0_rt);
+    vec_rubar[lev]->setVal(0.0_rt);
+    vec_rvbar[lev]->setVal(0.0_rt);
+    vec_rzeta[lev]->setVal(0.0_rt);
+
+    vec_ubar[lev]->setVal(0.0_rt);
+    vec_vbar[lev]->setVal(0.0_rt);
+    vec_zeta[lev]->setVal(0.0_rt);
 
     // Set initial linear drag coefficient
     vec_rdrag[lev]->setVal(solverChoice.rdrag);