diff --git a/Build/cmake.sh b/Build/cmake.sh index 4e49a5e5..5338ab09 100755 --- a/Build/cmake.sh +++ b/Build/cmake.sh @@ -13,6 +13,5 @@ cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \ -DROMSX_ENABLE_TESTS:BOOL=ON \ -DROMSX_ENABLE_FCOMPARE:BOOL=ON \ -DROMSX_ENABLE_DOCUMENTATION:BOOL=OFF \ - -DROMSX_ENABLE_SALINITY:BOOL=ON \ -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=ON \ .. && make -j8 diff --git a/CMake/BuildROMSXExe.cmake b/CMake/BuildROMSXExe.cmake index 9d36559e..8774af85 100644 --- a/CMake/BuildROMSXExe.cmake +++ b/CMake/BuildROMSXExe.cmake @@ -20,10 +20,6 @@ function(build_romsx_lib romsx_lib_name) ${ROMSX_EOS_DIR}/EOS.H) target_include_directories(${romsx_lib_name} SYSTEM PUBLIC ${ROMSX_EOS_DIR}) - if(ROMSX_ENABLE_SALINITY) - target_compile_definitions(${romsx_lib_name} PUBLIC ROMSX_USE_SALINITY) - endif() - if(ROMSX_ENABLE_PARTICLES) target_sources(${romsx_lib_name} PRIVATE ${SRC_DIR}/Particles/TracerPC.cpp) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8d478c2b..dc3c18d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +10,6 @@ include(CMakePackageConfigHelpers) #General options for all executables in the project set(ROMSX_DIM "3" CACHE STRING "Number of physical dimensions") option(ROMSX_ENABLE_DOCUMENTATION "Build documentation" OFF) -option(ROMSX_ENABLE_SALINITY "Enable Salinity" ON) option(ROMSX_ENABLE_ALL_WARNINGS "Enable all compiler warnings" OFF) option(ROMSX_ENABLE_TESTS "Enable regression and unit tests" OFF) option(ROMSX_ENABLE_NETCDF "Enable NetCDF IO" OFF) diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index 53d420f6..ed215978 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -2,4 +2,5 @@ set(romsx_lib_name romsx_srclib) add_library(${romsx_lib_name} OBJECT) include(${CMAKE_SOURCE_DIR}/CMake/BuildROMSXExe.cmake) build_romsx_lib(${romsx_lib_name}) +add_subdirectory(Seamount) add_subdirectory(Upwelling) diff --git a/Exec/Make.ROMSX b/Exec/Make.ROMSX index 24e3ab4e..46519494 100644 --- a/Exec/Make.ROMSX +++ b/Exec/Make.ROMSX @@ -73,6 +73,7 @@ endif #turn on NetCDF macro define ifeq ($(USE_NETCDF), TRUE) + USERSuffix = .NC DEFINES += -DROMSX_USE_NETCDF includes += $(shell pkg-config --cflags netcdf) LIBRARIES += $(shell pkg-config --libs netcdf) @@ -89,10 +90,6 @@ ifeq ($(USE_PARTICLES), TRUE) DEFINES += -DROMSX_USE_PARTICLES endif -ifeq ($(SALINITY),TRUE) - DEFINES += -DROMSX_USE_SALINITY -endif - CEXE_sources += AMReX_buildInfo.cpp CEXE_headers += $(AMREX_HOME)/Tools/C_scripts/AMReX_buildInfo.H INCLUDE_LOCATIONS += $(AMREX_HOME)/Tools/C_scripts diff --git a/Exec/ParticlesOverSeaMount/GNUmakefile b/Exec/ParticlesOverSeaMount/GNUmakefile index 655e7ff8..6bf59fec 100644 --- a/Exec/ParticlesOverSeaMount/GNUmakefile +++ b/Exec/ParticlesOverSeaMount/GNUmakefile @@ -26,9 +26,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -# Physics -SALINITY = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/ParticlesOverSeaMount/prob.cpp b/Exec/ParticlesOverSeaMount/prob.cpp index 9df31a3a..1ef68dac 100644 --- a/Exec/ParticlesOverSeaMount/prob.cpp +++ b/Exec/ParticlesOverSeaMount/prob.cpp @@ -72,7 +72,7 @@ init_custom_bathymetry (const Geometry& geom, Gpu::streamSynchronize(); - if(!m_solverChoice.flat_bathymetry) { + if (!m_solverChoice.flat_bathymetry) { ParallelFor(makeSlab(gbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { Real val1, val2; @@ -122,10 +122,12 @@ init_custom_prob( Array4 const& /*h*/, Array4 const& /*Zt_avg1*/, GeometryData const& geomdata, - const SolverChoice& /*m_solverChoice*/) + const SolverChoice& m_solverChoice) { const int khi = geomdata.Domain().bigEnd()[2]; + bool l_use_salt = m_solverChoice.use_salt; + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -141,9 +143,9 @@ init_custom_prob( state(i, j, k, Temp_comp) = 1.; state(i,j,k,Temp_comp)=parms.T0+8.0*std::exp(z/50.0_rt); -#ifdef ROMSX_USE_SALINITY - state(i,j,k,Salt_comp)=parms.S0; -#endif + if (l_use_salt) { + state(i,j,k,Salt_comp)=parms.S0; + } // Set scalar = 0 everywhere state(i, j, k, Scalar_comp) = 0.0; @@ -219,26 +221,27 @@ init_custom_vmix(const Geometry& /*geom*/, MultiFab& mf_Akv, MultiFab& mf_Akt, M void init_custom_hmix(const Geometry& /*geom*/, MultiFab& mf_visc2_p, MultiFab& mf_visc2_r, - MultiFab& mf_diff2_salt, MultiFab& mf_diff2_temp, const SolverChoice& /*m_solverChoice*/) + MultiFab& mf_diff2, const SolverChoice& /*m_solverChoice*/) { for ( MFIter mfi((mf_visc2_p), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& visc2_p = (mf_visc2_p).array(mfi); Array4 const& visc2_r = (mf_visc2_r).array(mfi); - Array4 const& diff2_salt = (mf_diff2_salt).array(mfi); - Array4 const& diff2_temp = (mf_diff2_temp).array(mfi); + Array4 const& diff2 = mf_diff2.array(mfi); Box bx = mfi.tilebox(); bx.grow(IntVect(NGROW,NGROW,0)); - Gpu::streamSynchronize(); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + int ncomp = mf_diff2.nComp(); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - visc2_p(i,j,k) = 5.0; - visc2_r(i,j,k) = 5.0; + visc2_p(i,j,k) = 5.0; + visc2_r(i,j,k) = 5.0; - diff2_salt(i,j,k) = 0.0; - diff2_temp(i,j,k) = 0.0; + for (int n = 0; n < ncomp; n++) { + diff2(i,j,k,n) = 0.0; + } }); } } diff --git a/Exec/Seamount/GNUmakefile b/Exec/Seamount/GNUmakefile index 0f01a5df..056c3f22 100644 --- a/Exec/Seamount/GNUmakefile +++ b/Exec/Seamount/GNUmakefile @@ -23,9 +23,6 @@ USE_ASSERTION = TRUE # Debugging DEBUG = TRUE -# Physics -SALINITY = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/Seamount/inputs b/Exec/Seamount/inputs index 64bab47e..c44b2213 100644 --- a/Exec/Seamount/inputs +++ b/Exec/Seamount/inputs @@ -49,7 +49,6 @@ romsx.plot_vars_1 = omega salt temp x_velocity y_velocity z_velocity romsx.plotfile_type = amrex # SOLVER CHOICE -romsx.use_gravity = true romsx.use_coriolis = true romsx.horizontal_advection_scheme = "upstream3" # upstream3 or centered4 romsx.spatial_order = 2 diff --git a/Exec/Seamount/prob.cpp b/Exec/Seamount/prob.cpp index 1e6c172e..7255802a 100644 --- a/Exec/Seamount/prob.cpp +++ b/Exec/Seamount/prob.cpp @@ -29,16 +29,7 @@ init_custom_bathymetry (const Geometry& geom, MultiFab& mf_Zt_avg1, const SolverChoice& m_solverChoice) { - //std::unique_ptr& mf_z_w = vec_z_w[lev]; - //std::unique_ptr& mf_h = vec_hOfTheConfusingName[lev]; - auto dx = geom.CellSizeArray(); - auto ProbLoArr = geom.ProbLoArray(); - auto ProbHiArr = geom.ProbHiArray(); - mf_h.setVal(geom.ProbHi(2)); - Real depth = geom.ProbHi(2); - const int Lm = geom.Domain().size()[0]; - const int Mm = geom.Domain().size()[1]; //HACK HACK manually setting zeta to 0 mf_Zt_avg1.setVal(0.0); @@ -54,50 +45,33 @@ init_custom_bathymetry (const Geometry& geom, const auto & geomdata = geom.data(); - int ncomp = 1; - //NOTE: this will need to be updated for multilevel - Box subdomain = geom.Domain(); - - int nx = subdomain.length(0); - int ny = subdomain.length(1); - int nz = subdomain.length(2); - - auto N = nz; // Number of vertical "levels" aka, NZ - bool NSPeriodic = geomdata.isPeriodic(1); - bool EWPeriodic = geomdata.isPeriodic(0); - amrex::Real Xsize = 320000.0_rt; amrex::Real Esize = 320000.0_rt; amrex::Real depth = 5000.0_rt; - if(!m_solverChoice.flat_bathymetry) { - Gpu::streamSynchronize(); - amrex::ParallelFor(gbx2, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + if(!m_solverChoice.flat_bathymetry) { - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); - - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real y = prob_lo[1] + (j + 0.5) * dx[1]; - - Real val1, val2; - int iFort = i+1; - int jFort = j+1; - val1 = (x-0.5_rt*Xsize)/40000.0_rt; - val2 = (y-0.5_rt*Esize)/40000.0_rt; - h(i,j,0) = depth - 4500.0_rt * std::exp(-(val1*val1+val2*val2)); - }); + ParallelFor(makeSlab(gbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int) + { + const auto prob_lo = geomdata.ProbLo(); + const auto dx = geomdata.CellSize(); + + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + + Real val1, val2; + val1 = (x-0.5_rt*Xsize)/40000.0_rt; + val2 = (y-0.5_rt*Esize)/40000.0_rt; + h(i,j,0) = depth - 4500.0_rt * std::exp(-(val1*val1+val2*val2)); + }); + } else { - Gpu::streamSynchronize(); - amrex::ParallelFor(gbx2, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) - { - h(i,j,0) = -geomdata.ProbLo(2); - if (k==0) { + + ParallelFor(makeSlab(gbx2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) + { + h(i,j,0) = -geomdata.ProbLo(2); h(i,j,0,1) = h(i,j,0); - } - }); + }); } } } @@ -111,38 +85,33 @@ init_custom_prob( Array4 const& z_vel, Array4 const& /*z_w*/, Array4 const& z_r, - Array4 const& Hz, - Array4 const& h, - Array4 const& Zt_avg1, + Array4 const& /*Hz*/, + Array4 const& /*h*/, + Array4 const& /*Zt_avg1*/, GeometryData const& geomdata, const SolverChoice& m_solverChoice) { - const int khi = geomdata.Domain().bigEnd()[2]; + bool l_use_salt = m_solverChoice.use_salt; - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + const int khi = geomdata.Domain().bigEnd()[2]; - const Real& rho_sfc = p_0 / (R_d*m_solverChoice.T0); - const Real& thetabar = m_solverChoice.T0; - const Real& dz = geomdata.CellSize()[2]; - const Real& el = geomdata.ProbHi()[1]; - const Real& prob_lo_z = geomdata.ProbLo()[2]; + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry (note we must include these here to get the data on device) - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); + // const auto prob_lo = geomdata.ProbLo(); + // const auto dx = geomdata.CellSize(); - const Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + // const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + // const Real y = prob_lo[1] + (j + 0.5) * dx[1]; const Real z = z_r(i,j,k); state(i, j, k, Temp_comp) = 1.; state(i,j,k,Temp_comp)=m_solverChoice.T0+7.5_rt*std::exp(z/1000.0_rt); -#ifdef ROMSX_USE_SALINITY - state(i,j,k,Salt_comp)=m_solverChoice.S0; -#endif + if (l_use_salt) { + state(i,j,k,Salt_comp)=m_solverChoice.S0; // Set scalar = 0 everywhere state(i, j, k, Scalar_comp) = 0.0; @@ -179,19 +148,20 @@ init_custom_prob( } void -init_custom_vmix(const Geometry& geom, MultiFab& mf_Akv, MultiFab& mf_Akt, MultiFab& mf_z_w, const SolverChoice& m_solverChoice) +init_custom_vmix(const Geometry& /*geom*/, MultiFab& mf_Akv, MultiFab& mf_Akt, + MultiFab& /*mf_z_w*/, const SolverChoice& /*m_solverChoice*/) { for ( MFIter mfi((mf_Akv), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& Akv = (mf_Akv).array(mfi); Array4 const& Akt = (mf_Akt).array(mfi); - Array4 const& z_w = (mf_z_w).array(mfi); + Box bx = mfi.tilebox(); bx.grow(IntVect(NGROW,NGROW,0)); - const auto & geomdata = geom.data(); - int ncomp = 1; Gpu::streamSynchronize(); - amrex::ParallelFor(bx, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + + Gpu::streamSynchronize(); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Akv(i,j,k) = 1.0e-5_rt; Akt(i,j,k) = 1.0e-6_rt; @@ -200,35 +170,36 @@ init_custom_vmix(const Geometry& geom, MultiFab& mf_Akv, MultiFab& mf_Akt, Multi } void -init_custom_hmix(const Geometry& geom, MultiFab& mf_visc2_p, MultiFab& mf_visc2_r, - MultiFab& mf_diff2_salt, MultiFab& mf_diff2_temp, const SolverChoice& m_solverChoice) +init_custom_hmix(const Geometry& /*geom*/, MultiFab& mf_visc2_p, MultiFab& mf_visc2_r, + MultiFab& mf_diff2, const SolverChoice& /*m_solverChoice*/) { for ( MFIter mfi((mf_visc2_p), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - Array4 const& visc2_p = (mf_visc2_p).array(mfi); - Array4 const& visc2_r = (mf_visc2_r).array(mfi); - Array4 const& diff2_salt = (mf_diff2_salt).array(mfi); - Array4 const& diff2_temp = (mf_diff2_temp).array(mfi); + Array4 const& visc2_p_arr = mf_visc2_p.array(mfi); + Array4 const& visc2_r_arr = mf_visc2_r.array(mfi); + Array4 const& diff2_arr = mf_diff2.array(mfi); + Box bx = mfi.tilebox(); bx.grow(IntVect(NGROW,NGROW,0)); - const auto & geomdata = geom.data(); - int ncomp = 1; + + int ncomp = 2; // temperature and salt Gpu::streamSynchronize(); - amrex::ParallelFor(bx, ncomp, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - visc2_p(i,j,k) = 0.0; - visc2_r(i,j,k) = 0.0; + visc2_p_arr(i,j,k) = 0.0; + visc2_r_arr(i,j,k) = 0.0; - diff2_salt(i,j,k) = 0.0; - diff2_temp(i,j,k) = 0.0; + for (int n = 0; n < ncomp; n++) { + diff2_arr(i,j,k,n) = 0.0; + } }); - } + } // mfi } void -init_custom_smflux(const Geometry& geom, const Real time, MultiFab& mf_sustr, MultiFab& mf_svstr, - const SolverChoice& m_solverChoice) +init_custom_smflux(const Geometry& /*geom*/, const Real /*time*/, + MultiFab& mf_sustr, MultiFab& mf_svstr, const SolverChoice& /*m_solverChoice*/) { mf_sustr.setVal(0.0); mf_svstr.setVal(0.0); diff --git a/Exec/Upwelling/GNUmakefile b/Exec/Upwelling/GNUmakefile index 6274bde3..6fc8716a 100644 --- a/Exec/Upwelling/GNUmakefile +++ b/Exec/Upwelling/GNUmakefile @@ -26,9 +26,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -# Physics -SALINITY = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/Upwelling/prob.cpp b/Exec/Upwelling/prob.cpp index a893f2ef..85f32e4e 100644 --- a/Exec/Upwelling/prob.cpp +++ b/Exec/Upwelling/prob.cpp @@ -119,9 +119,11 @@ init_custom_prob( GeometryData const& geomdata, const SolverChoice& m_solverChoice) { - const int khi = geomdata.Domain().bigEnd()[2]; + bool l_use_salt = m_solverChoice.use_salt; - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); + const int khi = geomdata.Domain().bigEnd()[2]; + + AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -136,9 +138,9 @@ init_custom_prob( state(i, j, k, Temp_comp) = 1.; state(i,j,k,Temp_comp)=m_solverChoice.T0+8.0*std::exp(z/50.0_rt); -#ifdef ROMSX_USE_SALINITY - state(i,j,k,Salt_comp)=m_solverChoice.S0; -#endif + if (l_use_salt) { + state(i,j,k,Salt_comp)=m_solverChoice.S0; + } // Set scalar = 0 everywhere state(i, j, k, Scalar_comp) = 0.0; @@ -211,25 +213,27 @@ init_custom_vmix(const Geometry& /*geom*/, MultiFab& mf_Akv, MultiFab& mf_Akt, void init_custom_hmix(const Geometry& /*geom*/, MultiFab& mf_visc2_p, MultiFab& mf_visc2_r, - MultiFab& mf_diff2_salt, MultiFab& mf_diff2_temp, const SolverChoice& /*m_solverChoice*/) + MultiFab& mf_diff2, const SolverChoice& /*m_solverChoice*/) { for ( MFIter mfi((mf_visc2_p), TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& visc2_p = (mf_visc2_p).array(mfi); Array4 const& visc2_r = (mf_visc2_r).array(mfi); - Array4 const& diff2_salt = (mf_diff2_salt).array(mfi); - Array4 const& diff2_temp = (mf_diff2_temp).array(mfi); + Array4 const& diff2 = mf_diff2.array(mfi); Box bx = mfi.tilebox(); bx.grow(IntVect(NGROW,NGROW,0)); Gpu::streamSynchronize(); - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + + int ncomp = mf_diff2.nComp(); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { visc2_p(i,j,k) = 5.0; visc2_r(i,j,k) = 5.0; - diff2_salt(i,j,k) = 0.0; - diff2_temp(i,j,k) = 0.0; + for (int n = 0; n < ncomp; n++) { + diff2(i,j,k,n) = 0.0; + } }); } } diff --git a/Source/BoundaryConditions/BoundaryConditions_cons.cpp b/Source/BoundaryConditions/BoundaryConditions_cons.cpp index cea342f7..324d1b59 100644 --- a/Source/BoundaryConditions/BoundaryConditions_cons.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_cons.cpp @@ -6,16 +6,16 @@ using namespace amrex; // // mf is the multifab to be filled // icomp is the index into the MultiFab -- if cell-centered this can be any value -// from 0 to NVAR-1, if face-centered this must be 0 +// from 0 to NCONS-1, if face-centered this must be 0 // ncomp is the number of components -- if cell-centered (var_idx = 0) this can be any value -// from 1 to NVAR as long as icomp+ncomp <= NVAR-1. If face-centered this +// from 1 to NCONS as long as icomp+ncomp <= NCONS-1. If face-centered this // must be 1 // time is the time at which the data should be filled // bccomp is the index into both domain_bcs_type_bcr and bc_extdir_vals for icomp = 0 -- // so this follows the BCVars enum // void ROMSXPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box& bx, const Box& domain, - const GpuArray dxInv, + const GpuArray /*dxInv*/, int icomp, int ncomp, Real /*time*/, int bccomp) { BL_PROFILE_VAR("impose_cons_bcs()",impose_cons_bcs); @@ -43,7 +43,7 @@ void ROMSXPhysBCFunct::impose_cons_bcs (const Array4& dest_arr, const Box& #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NCONS> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp index 3320f17c..4dcea676 100644 --- a/Source/BoundaryConditions/BoundaryConditions_xvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_xvel.cpp @@ -39,7 +39,7 @@ void ROMSXPhysBCFunct::impose_xvel_bcs (const Array4& dest_arr, const Box& #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NCONS> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp index b38758a2..b0af20e2 100644 --- a/Source/BoundaryConditions/BoundaryConditions_yvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_yvel.cpp @@ -39,7 +39,7 @@ void ROMSXPhysBCFunct::impose_yvel_bcs (const Array4& dest_arr, const Box& #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray, AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray, AMREX_SPACEDIM+NCONS> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp index 2536f259..7ed3a7e1 100644 --- a/Source/BoundaryConditions/BoundaryConditions_zvel.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_zvel.cpp @@ -47,7 +47,7 @@ void ROMSXPhysBCFunct::impose_zvel_bcs (const Array4& dest_arr, const Box& #endif const amrex::BCRec* bc_ptr = bcrs_d.data(); - GpuArray,AMREX_SPACEDIM+NVAR> l_bc_extdir_vals_d; + GpuArray,AMREX_SPACEDIM+NCONS> l_bc_extdir_vals_d; for (int i = 0; i < ncomp; i++) for (int ori = 0; ori < 2*AMREX_SPACEDIM; ori++) diff --git a/Source/BoundaryConditions/ROMSX_FillPatch.cpp b/Source/BoundaryConditions/ROMSX_FillPatch.cpp index 68ea41ba..25ada60d 100644 --- a/Source/BoundaryConditions/ROMSX_FillPatch.cpp +++ b/Source/BoundaryConditions/ROMSX_FillPatch.cpp @@ -166,7 +166,7 @@ ROMSX::FillCoarsePatch (int lev, Real time, MultiFab* mf_to_fill, MultiFab* mf_c { bccomp = 0; mapper = &cell_cons_interp; - ncomp = Cons::NumVars; + ncomp = NCONS; } else if (box_mf.ixType() == IndexType(IntVect(1,0,0))) { diff --git a/Source/BoundaryConditions/ROMSX_PhysBCFunct.H b/Source/BoundaryConditions/ROMSX_PhysBCFunct.H index 66d0c73a..f711ded8 100644 --- a/Source/BoundaryConditions/ROMSX_PhysBCFunct.H +++ b/Source/BoundaryConditions/ROMSX_PhysBCFunct.H @@ -37,7 +37,7 @@ public: const amrex::Geometry& geom, const amrex::Vector& domain_bcs_type, const amrex::Gpu::DeviceVector& domain_bcs_type_d, TimeInterpolatedData& data, - amrex::Array,AMREX_SPACEDIM+NVAR> bc_extdir_vals + amrex::Array,AMREX_SPACEDIM+NCONS> bc_extdir_vals #ifdef ROMSX_USE_NETCDF ,const std::string& init_type, const amrex::Vector>& bdy_data_xlo, @@ -66,9 +66,9 @@ public: // // mf is the multifab to be filled // icomp is the index into the MultiFab -- if cell-centered this can be any value - // from 0 to NVAR-1, if face-centered this must be 0 + // from 0 to NCONS-1, if face-centered this must be 0 // ncomp is the number of components -- if cell-centered this can be any value - // from 1 to NVAR as long as icomp+ncomp <= NVAR-1. If face-centered this + // from 1 to NCONS as long as icomp+ncomp <= NCONS-1. If face-centered this // must be 1 // nghost is how many ghost cells to be filled // time is the time at which the data should be filled @@ -111,7 +111,7 @@ private: amrex::Vector m_domain_bcs_type; amrex::Gpu::DeviceVector m_domain_bcs_type_d; TimeInterpolatedData& m_data; - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NCONS> m_bc_extdir_vals; #ifdef ROMSX_USE_NETCDF const std::string m_init_type; const amrex::Vector>& m_bdy_data_xlo; diff --git a/Source/BoundaryConditions/ROMSX_PhysBCFunct.cpp b/Source/BoundaryConditions/ROMSX_PhysBCFunct.cpp index af678c58..e871b9be 100644 --- a/Source/BoundaryConditions/ROMSX_PhysBCFunct.cpp +++ b/Source/BoundaryConditions/ROMSX_PhysBCFunct.cpp @@ -7,9 +7,9 @@ using namespace amrex; // // mf is the multifab to be filled // icomp is the index into the MultiFab -- if cell-centered this can be any value -// from 0 to NVAR-1, if face-centered this must be 0 +// from 0 to NCONS-1, if face-centered this must be 0 // ncomp is the number of components -- if cell-centered this can be any value -// from 1 to NVAR as long as icomp+ncomp <= NVAR-1. If face-centered this +// from 1 to NCONS as long as icomp+ncomp <= NCONS-1. If face-centered this // must be 1 // nghost is how many ghost cells to be filled // time is the time at which the data should be filled @@ -76,7 +76,7 @@ void ROMSXPhysBCFunct::operator() (MultiFab& mf, int icomp, int ncomp, IntVect c } else if (mf[0].box().ixType() == IndexType(IntVect(0,0,0))) { - AMREX_ALWAYS_ASSERT(icomp == 0 && icomp+ncomp <= NVAR); + AMREX_ALWAYS_ASSERT(icomp == 0 && icomp+ncomp <= NCONS); impose_cons_bcs(dest_arr,bx,domain, dxInv,icomp,ncomp,time,bccomp); } else { diff --git a/Source/DataStruct.H b/Source/DataStruct.H index 4eb9fbdd..20799eaa 100644 --- a/Source/DataStruct.H +++ b/Source/DataStruct.H @@ -53,6 +53,9 @@ struct SolverChoice { // Order of spatial discretization pp.query("spatial_order", spatial_order); + // Include salinity? + pp.query("use_salt", use_salt); + // Include Coriolis forcing? pp.query("use_coriolis", use_coriolis); @@ -111,6 +114,7 @@ struct SolverChoice { void display() { amrex::Print() << "SOLVER CHOICE: " << std::endl; + amrex::Print() << "use_salt : " << use_salt << std::endl; amrex::Print() << "use_coriolis : " << use_coriolis << std::endl; amrex::Print() << "use_prestep : " << use_prestep << std::endl; amrex::Print() << "use_uv3dmix : " << use_uv3dmix << std::endl; @@ -147,6 +151,8 @@ struct SolverChoice { bool flat_bathymetry = true; + bool use_salt = true; + // Specify what additional physics/forcing modules we use bool use_coriolis = false; diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp index cfebdb59..7c897632 100644 --- a/Source/IO/Checkpoint.cpp +++ b/Source/IO/Checkpoint.cpp @@ -62,7 +62,7 @@ ROMSX::WriteCheckpointFile () const // for each variable we store // conservative, cell-centered vars - HeaderFile << Cons::NumVars << "\n"; + HeaderFile << NCONS << "\n"; // x-velocity on faces HeaderFile << 1 << "\n"; @@ -116,8 +116,8 @@ ROMSX::WriteCheckpointFile () const } BoxArray ba2d(std::move(bl2d)); - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); - MultiFab::Copy(cons,*cons_new[lev],0,0,NVAR,0); + MultiFab cons(grids[lev],dmap[lev],NCONS,0); + MultiFab::Copy(cons,*cons_new[lev],0,0,NCONS,0); VisMF::Write(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell")); MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); @@ -205,7 +205,7 @@ ROMSX::ReadCheckpointFile () // conservative, cell-centered vars is >> chk_ncomp; GotoNextLine(is); - AMREX_ASSERT(chk_ncomp == Cons::NumVars); + AMREX_ASSERT(chk_ncomp == NCONS); // x-velocity on faces is >> chk_ncomp; @@ -290,9 +290,9 @@ ROMSX::ReadCheckpointFile () } BoxArray ba2d(std::move(bl2d)); - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); + MultiFab cons(grids[lev],dmap[lev],NCONS,0); VisMF::Read(cons, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell")); - MultiFab::Copy(*cons_new[lev],cons,0,0,NVAR,0); + MultiFab::Copy(*cons_new[lev],cons,0,0,NCONS,0); MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); VisMF::Read(xvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "XFace")); diff --git a/Source/IO/NCCheckpoint.cpp b/Source/IO/NCCheckpoint.cpp index 803f6c19..9e0fd84b 100644 --- a/Source/IO/NCCheckpoint.cpp +++ b/Source/IO/NCCheckpoint.cpp @@ -97,8 +97,8 @@ ROMSX::WriteNCCheckpointFile () const // Here we make copies of the MultiFab with no ghost cells for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); - MultiFab::Copy(cons,*cons_new[lev],0,0,NVAR,0); + MultiFab cons(grids[lev],dmap[lev],NCONS,0); + MultiFab::Copy(cons,*cons_new[lev],0,0,NCONS,0); WriteNCMultiFab(cons, amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Cell")); MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); @@ -136,7 +136,7 @@ ROMSX::ReadNCCheckpointFile () const std::string ntime_name = "num_newtime"; const int nvar = static_cast(ncf.dim(nvar_name).len()); - AMREX_ALWAYS_ASSERT(nvar == Cons::NumVars); + AMREX_ALWAYS_ASSERT(nvar == NCONS); const int ndt = static_cast(ncf.dim(ndt_name).len()); const int nstep = static_cast(ncf.dim(nstep_name).len()); @@ -188,7 +188,7 @@ ROMSX::ReadNCCheckpointFile () SetDistributionMap(lev, dm); // build MultiFab data - int ncomp = Cons::NumVars; + int ncomp = NCONS; cons_new[lev]->define(grids[lev], dmap[lev], ncomp, ngrow_state); cons_old[lev]->define(grids[lev], dmap[lev], ncomp, ngrow_state); @@ -207,9 +207,9 @@ ROMSX::ReadNCCheckpointFile () for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab cons(grids[lev],dmap[lev],Cons::NumVars,0); + MultiFab cons(grids[lev],dmap[lev],NCONS,0); WriteNCMultiFab(cons, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell")); - MultiFab::Copy(*cons_new[lev],cons,0,0,Cons::NumVars,0); + MultiFab::Copy(*cons_new[lev],cons,0,0,NCONS,0); MultiFab xvel(convert(grids[lev],IntVect(1,0,0)),dmap[lev],1,0); WriteNCMultiFab(xvel, amrex::MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Cell")); @@ -224,7 +224,7 @@ ROMSX::ReadNCCheckpointFile () MultiFab::Copy(*zvel_new[lev],zvel,0,0,1,0); // Copy from new into old just in case - MultiFab::Copy(*cons_old[lev],*cons_new[lev],0,0,NVAR,0); + MultiFab::Copy(*cons_old[lev],*cons_new[lev],0,0,NCONS,0); MultiFab::Copy(*xvel_old[lev],*xvel_new[lev],0,0,1,0); MultiFab::Copy(*yvel_old[lev],*yvel_new[lev],0,0,1,0); MultiFab::Copy(*zvel_old[lev],*zvel_new[lev],0,0,1,0); diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index d3f0330e..ddf85acf 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -37,7 +37,7 @@ ROMSX::setPlotVariables (const std::string& pp_plot_var_names, Vector tmp_plot_names; - for (int i = 0; i < Cons::NumVars; ++i) { + for (int i = 0; i < NCONS; ++i) { if ( containerHasElement(plot_var_names, cons_names[i]) ) { tmp_plot_names.push_back(cons_names[i]); } @@ -115,8 +115,8 @@ ROMSX::WritePlotFile (int which, Vector plot_var_names) int mf_comp = 0; // First, copy any of the conserved state variables into the output plotfile - AMREX_ALWAYS_ASSERT(cons_names.size() == Cons::NumVars); - for (int i = 0; i < Cons::NumVars; ++i) { + AMREX_ALWAYS_ASSERT(cons_names.size() == NCONS); + for (int i = 0; i < NCONS; ++i) { if (containerHasElement(plot_var_names, cons_names[i])) { MultiFab::Copy(mf[lev],*cons_new[lev],i,mf_comp,1,ngrow_vars); mf_comp++; diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 8efa17e7..44cbb441 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -8,13 +8,8 @@ #define Temp_comp 0 #define Scalar_comp 1 #define Omega_comp 2 - -#ifdef ROMSX_USE_SALINITY #define Salt_comp 3 -#define NVAR 4 -#else -#define NVAR 3 -#endif +#define NCONS 4 #define NGROW 2 @@ -23,25 +18,13 @@ namespace BCVars { cons_bc = 0, Temp_bc_comp = 0, Scalar_bc_comp, - xvel_bc = NVAR, + xvel_bc = NCONS, yvel_bc, zvel_bc, NumTypes }; } -namespace Cons { - enum { - Temp = 0, - Scalar, - Omega, -#ifdef ROMSX_USE_SALINITY - Salt, -#endif - NumVars - }; -} - enum struct ROMSX_BC { symmetry, inflow, outflow, no_slip_wall, slip_wall, periodic, MOST, undefined }; diff --git a/Source/Initialization/ROMSX_init_bcs.cpp b/Source/Initialization/ROMSX_init_bcs.cpp index fd40c331..9835b386 100644 --- a/Source/Initialization/ROMSX_init_bcs.cpp +++ b/Source/Initialization/ROMSX_init_bcs.cpp @@ -119,8 +119,8 @@ void ROMSX::init_bcs () // // ***************************************************************************** { - domain_bcs_type.resize(AMREX_SPACEDIM+NVAR); - domain_bcs_type_d.resize(AMREX_SPACEDIM+NVAR); + domain_bcs_type.resize(AMREX_SPACEDIM+NCONS); + domain_bcs_type_d.resize(AMREX_SPACEDIM+NCONS); for (OrientationIter oit; oit; ++oit) { Orientation ori = oit(); @@ -214,32 +214,32 @@ void ROMSX::init_bcs () if ( bct == ROMSX_BC::symmetry ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ROMSXBCType::reflect_even); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ROMSXBCType::reflect_even); } } else if ( bct == ROMSX_BC::outflow ) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ROMSXBCType::foextrap); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ROMSXBCType::foextrap); } } else if ( bct == ROMSX_BC::no_slip_wall) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ROMSXBCType::foextrap); if (m_bc_extdir_vals[BCVars::Temp_bc_comp][ori] > 0.) domain_bcs_type[BCVars::Temp_bc_comp].setLo(dir, ROMSXBCType::ext_dir); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ROMSXBCType::foextrap); if (m_bc_extdir_vals[BCVars::Temp_bc_comp][ori] > 0.) domain_bcs_type[BCVars::Temp_bc_comp].setHi(dir, ROMSXBCType::ext_dir); @@ -248,12 +248,12 @@ void ROMSX::init_bcs () else if (bct == ROMSX_BC::slip_wall) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ROMSXBCType::foextrap); if (m_bc_extdir_vals[BCVars::Temp_bc_comp][ori] > 0.) domain_bcs_type[BCVars::Temp_bc_comp].setLo(dir, ROMSXBCType::ext_dir); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ROMSXBCType::foextrap); if (m_bc_extdir_vals[BCVars::Temp_bc_comp][ori] > 0.) domain_bcs_type[BCVars::Temp_bc_comp].setHi(dir, ROMSXBCType::ext_dir); @@ -262,11 +262,11 @@ void ROMSX::init_bcs () else if (bct == ROMSX_BC::inflow) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) { + for (int i = 0; i < NCONS; i++) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ROMSXBCType::ext_dir); } } else { - for (int i = 0; i < NVAR; i++) { + for (int i = 0; i < NCONS; i++) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ROMSXBCType::ext_dir); } } @@ -274,10 +274,10 @@ void ROMSX::init_bcs () else if (bct == ROMSX_BC::periodic) { if (side == Orientation::low) { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ROMSXBCType::int_dir); } else { - for (int i = 0; i < NVAR; i++) + for (int i = 0; i < NCONS; i++) domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ROMSXBCType::int_dir); } } @@ -287,11 +287,11 @@ void ROMSX::init_bcs () #ifdef AMREX_USE_GPU Gpu::htod_memcpy (domain_bcs_type_d.data(), domain_bcs_type.data(), - sizeof(amrex::BCRec)*(NVAR+AMREX_SPACEDIM)); + sizeof(amrex::BCRec)*(NCONS+AMREX_SPACEDIM)); #else std::memcpy (domain_bcs_type_d.data(), domain_bcs_type.data(), - sizeof(amrex::BCRec)*(NVAR+AMREX_SPACEDIM)); + sizeof(amrex::BCRec)*(NCONS+AMREX_SPACEDIM)); #endif } diff --git a/Source/Initialization/prob_common.H b/Source/Initialization/prob_common.H index e8412353..19f384fe 100644 --- a/Source/Initialization/prob_common.H +++ b/Source/Initialization/prob_common.H @@ -39,8 +39,7 @@ void init_custom_vmix (const amrex::Geometry& geom, void init_custom_hmix (const amrex::Geometry& geom, amrex::MultiFab& mf_visc2_p, amrex::MultiFab& mf_visc2_r, - amrex::MultiFab& mf_diff2_salt, - amrex::MultiFab& mf_diff2_temp, + amrex::MultiFab& mf_diff2, SolverChoice const& m_solverChoice); void init_custom_smflux (const amrex::Geometry& geom, diff --git a/Source/Particles/TracerPC.cpp b/Source/Particles/TracerPC.cpp index 187ad9f4..3287af9d 100644 --- a/Source/Particles/TracerPC.cpp +++ b/Source/Particles/TracerPC.cpp @@ -176,14 +176,14 @@ TracerPC::AdvectWithUmac (Array umac, &((*umac_pointer[2])[grid])) }; //array of these pointers to pass to the GPU - amrex::GpuArray, AMREX_SPACEDIM> + GpuArray, AMREX_SPACEDIM> const umacarr {{AMREX_D_DECL((*fab[0]).array(), (*fab[1]).array(), (*fab[2]).array() )}}; auto zheight = use_terrain ? a_z_height[grid].array() : Array4{}; - amrex::ParallelFor(n, + ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; diff --git a/Source/ROMSX.H b/Source/ROMSX.H index 3ec17d0d..b2ec3665 100644 --- a/Source/ROMSX.H +++ b/Source/ROMSX.H @@ -215,10 +215,8 @@ public: amrex::Vector> vec_visc2_p; /** Harmonic viscosity devined on the rho points (centers) */ amrex::Vector> vec_visc2_r; - /** Harmonic diffusivity for salinity */ - amrex::Vector> vec_diff2_salt; - /** Harmonic diffusivity for temperature */ - amrex::Vector> vec_diff2_temp; + /** Harmonic diffusivity for temperature / salinity */ + amrex::Vector> vec_diff2; /** x coordinates at rho points (cell centers) */ amrex::Vector> vec_x_r; @@ -314,7 +312,8 @@ public: std::unique_ptr& mf_h, std::unique_ptr& mf_visc2_p, std::unique_ptr& mf_visc2_r, - const int ncomp, const amrex::Real dtfast_lev, + const int ncomp, + const amrex::Real dtfast_lev, const bool predictor_2d_step, const bool first_2d_step, int my_iif, int & next_indx1); @@ -338,7 +337,7 @@ public: std::unique_ptr& mf_Hvom, std::unique_ptr& mf_z_w, std::unique_ptr& mf_h , - const int ncomp, const int N, + const int N, const amrex::Real dt_lev); /** Wrapper function for prestep */ @@ -347,8 +346,8 @@ public: amrex::MultiFab& mf_u, amrex::MultiFab& mf_v, std::unique_ptr& mf_ru, std::unique_ptr& mf_rv, - amrex::MultiFab& mf_tempold, amrex::MultiFab& mf_saltold, - amrex::MultiFab& mf_temp, amrex::MultiFab& mf_salt, + amrex::MultiFab& S_old, + amrex::MultiFab& S_new, amrex::MultiFab& mf_W, amrex::MultiFab& mf_DC, /* MF mf_FC? */ std::unique_ptr& mf_z_r, @@ -396,10 +395,6 @@ public: amrex::Array4 rv, amrex::Array4 Hz, amrex::Array4 Akv, - amrex::Array4 on_u, - amrex::Array4 om_v, - amrex::Array4 Huon, - amrex::Array4 Hvom, amrex::Array4 pm, amrex::Array4 pn, amrex::Array4 W, @@ -461,15 +456,14 @@ public: int nrhs); void rho_eos (const amrex::Box& bx, - amrex::Array4 temp, - amrex::Array4 salt, + amrex::Array4 state, amrex::Array4 rho, amrex::Array4 rhoA, amrex::Array4 rhoS, amrex::Array4 Hz, amrex::Array4 z_w, amrex::Array4 h, - const int nrhs, const int N); + const int N); void prsgrd (const amrex::Box& bx, const amrex::Box& gbx, @@ -526,8 +520,7 @@ public: amrex::Array4 v, amrex::Array4 Hv, amrex::Array4 om_v, - amrex::Array4 Hz, - const int nnew); + amrex::Array4 Hz); void update_massflux_3d (const amrex::Box& bx, const int ioff, const int joff, @@ -575,15 +568,15 @@ public: /** Harmonic diffusivity for tracers */ void t3dmix (const amrex::Box& bx, - amrex::Array4 t, + amrex::Array4 state, amrex::Array4 diff2, amrex::Array4 Hz, amrex::Array4 pm, amrex::Array4 pn, amrex::Array4 pmon_u, amrex::Array4 pnom_v, - int nrhs, int nnew, - const amrex::Real dt_lev); + const amrex::Real dt_lev, + const int ncomp); void coriolis (const amrex::Box& xbx, const amrex::Box& ybx, @@ -784,7 +777,7 @@ private: amrex::Array domain_bc_type; // These hold the Dirichlet values at walls which need them ... - amrex::Array,AMREX_SPACEDIM+NVAR> m_bc_extdir_vals; + amrex::Array,AMREX_SPACEDIM+NCONS> m_bc_extdir_vals; // These are the "physical" boundary condition types (e.g. "inflow") amrex::GpuArray phys_bc_type; diff --git a/Source/ROMSX.cpp b/Source/ROMSX.cpp index f95dac1d..f2f70387 100644 --- a/Source/ROMSX.cpp +++ b/Source/ROMSX.cpp @@ -49,7 +49,7 @@ amrex::Vector BCNames = {"xlo", "ylo", "zlo", "xhi", "yhi", "zhi"}; // - initializes BCRec boundary condition object ROMSX::ROMSX () { - if (amrex::ParallelDescriptor::IOProcessor()) { + if (ParallelDescriptor::IOProcessor()) { const char* romsx_hash = amrex::buildInfoGetGitHash(1); const char* amrex_hash = amrex::buildInfoGetGitHash(2); const char* buildgithash = amrex::buildInfoGetBuildGitHash(); @@ -222,7 +222,7 @@ ROMSX::post_timestep (int nstep, Real time, Real dt_lev0) for (int lev = finest_level-1; lev >= 0; lev--) { // This call refluxes from the lev/lev+1 interface onto lev - getAdvFluxReg(lev+1)->Reflux(*cons_new[lev], 0, 0, NVAR); + getAdvFluxReg(lev+1)->Reflux(*cons_new[lev], 0, 0, NCONS); // We need to do this before anything else because refluxing changes the // values of coarse cells underneath fine grids with the assumption they'll @@ -241,7 +241,7 @@ void ROMSX::InitData () { // Initialize the start time for our CPU-time tracker - startCPUTime = amrex::ParallelDescriptor::second(); + startCPUTime = ParallelDescriptor::second(); // Map the words in the inputs file to BC types, then translate // those types into what they mean for each variable @@ -265,7 +265,7 @@ ROMSX::InitData () } #ifdef ROMSX_USE_PARTICLES - particleData.init_particles((amrex::ParGDBBase*)GetParGDB(),vec_z_phys_nd); + particleData.init_particles((ParGDBBase*)GetParGDB(),vec_z_phys_nd); #endif } else { // Restart from a checkpoint @@ -282,7 +282,7 @@ ROMSX::InitData () advflux_reg[lev] = new YAFluxRegister(grids[lev], grids[lev-1], dmap[lev], dmap[lev-1], geom[lev], geom[lev-1], - ref_ratio[lev-1], lev, NVAR); + ref_ratio[lev-1], lev, NCONS); } } @@ -333,7 +333,7 @@ ROMSX::InitData () // int ngs = cons_new[lev]->nGrow(); int ngvel = xvel_new[lev]->nGrow(); - MultiFab::Copy(*cons_old[lev],*cons_new[lev],0,0,NVAR,ngs); + MultiFab::Copy(*cons_old[lev],*cons_new[lev],0,0,NCONS,ngs); MultiFab::Copy(*xvel_old[lev],*xvel_new[lev],0,0,1,ngvel); MultiFab::Copy(*yvel_old[lev],*yvel_new[lev],0,0,1,ngvel); MultiFab::Copy(*zvel_old[lev],*zvel_new[lev],0,0,1,IntVect(ngvel,ngvel,0)); @@ -375,13 +375,11 @@ ROMSX::set_vmix(int lev) { void ROMSX::set_hmixcoef(int lev) { - init_custom_hmix(geom[lev], *vec_visc2_p[lev], *vec_visc2_r[lev], *vec_diff2_salt[lev], - *vec_diff2_temp[lev], solverChoice); + init_custom_hmix(geom[lev], *vec_visc2_p[lev], *vec_visc2_r[lev], *vec_diff2[lev], solverChoice); vec_visc2_p[lev]->FillBoundary(geom[lev].periodicity()); vec_visc2_r[lev]->FillBoundary(geom[lev].periodicity()); - vec_diff2_salt[lev]->FillBoundary(geom[lev].periodicity()); - vec_diff2_temp[lev]->FillBoundary(geom[lev].periodicity()); + vec_diff2[lev]->FillBoundary(geom[lev].periodicity()); } void diff --git a/Source/ROMSX_SumIQ.cpp b/Source/ROMSX_SumIQ.cpp index 5aabbdc8..c3eab4e4 100644 --- a/Source/ROMSX_SumIQ.cpp +++ b/Source/ROMSX_SumIQ.cpp @@ -16,8 +16,8 @@ ROMSX::sum_integrated_quantities(Real time) int datwidth = 14; int datprecision = 6; - amrex::Real scalar = 0.0; - amrex::Real kineng = 0.0; + Real scalar = 0.0; + Real kineng = 0.0; for (int lev = 0; lev <= finest_level; lev++) { @@ -43,14 +43,14 @@ ROMSX::sum_integrated_quantities(Real time) if (verbose > 0) { const int nfoo = 2; - amrex::Real foo[nfoo] = {scalar,kineng}; + Real foo[nfoo] = {scalar,kineng}; #ifdef AMREX_LAZY Lazy::QueueReduction([=]() mutable { #endif - amrex::ParallelDescriptor::ReduceRealSum( - foo, nfoo, amrex::ParallelDescriptor::IOProcessorNumber()); + ParallelDescriptor::ReduceRealSum( + foo, nfoo, ParallelDescriptor::IOProcessorNumber()); - if (amrex::ParallelDescriptor::IOProcessor()) { + if (ParallelDescriptor::IOProcessor()) { int i = 0; scalar = foo[i++]; kineng = foo[i++]; @@ -86,8 +86,7 @@ ROMSX::sum_integrated_quantities(Real time) } Real -ROMSX::volWgtSumMF(int lev, - const amrex::MultiFab& mf, int comp, bool local, bool finemask) +ROMSX::volWgtSumMF(int lev, const MultiFab& mf, int comp, bool local, bool finemask) { BL_PROFILE("ROMSX::volWgtSumMF()"); @@ -104,7 +103,7 @@ ROMSX::volWgtSumMF(int lev, auto const& dx = geom[lev].CellSizeArray(); Real cell_vol = dx[0]*dx[1]*dx[2]; volume.setVal(cell_vol); - sum = amrex::MultiFab::Dot(tmp, 0, volume, 0, 1, 0, local); + sum = MultiFab::Dot(tmp, 0, volume, 0, 1, 0, local); if (!local) ParallelDescriptor::ReduceRealSum(sum); @@ -112,7 +111,7 @@ ROMSX::volWgtSumMF(int lev, return sum; } -amrex::MultiFab& +MultiFab& ROMSX::build_fine_mask(int level) { // Mask for zeroing covered cells @@ -122,26 +121,26 @@ ROMSX::build_fine_mask(int level) const DistributionMapping& cdm = dmap[level-1]; // TODO -- we should make a vector of these a member of ROMSX class - fine_mask.define(cba, cdm, 1, 0, amrex::MFInfo()); + fine_mask.define(cba, cdm, 1, 0, MFInfo()); fine_mask.setVal(1.0); - amrex::BoxArray fba = grids[level]; - amrex::iMultiFab ifine_mask = makeFineMask(cba, cdm, fba, ref_ratio[level-1], 1, 0); + BoxArray fba = grids[level]; + iMultiFab ifine_mask = makeFineMask(cba, cdm, fba, ref_ratio[level-1], 1, 0); const auto fma = fine_mask.arrays(); const auto ifma = ifine_mask.arrays(); - amrex::ParallelFor(fine_mask, [=] AMREX_GPU_DEVICE(int bno, int i, int j, int k) noexcept + ParallelFor(fine_mask, [=] AMREX_GPU_DEVICE(int bno, int i, int j, int k) noexcept { fma[bno](i,j,k) = ifma[bno](i,j,k); }); - amrex::Gpu::synchronize(); + Gpu::synchronize(); return fine_mask; } bool -ROMSX::is_it_time_for_action(int nstep, Real time, Real dtlev, int action_interval, amrex::Real action_per) +ROMSX::is_it_time_for_action(int nstep, Real time, Real dtlev, int action_interval, Real action_per) { bool int_test = (action_interval > 0 && nstep % action_interval == 0); diff --git a/Source/ROMSX_make_new_level.cpp b/Source/ROMSX_make_new_level.cpp index c3203301..1bbc3e0d 100644 --- a/Source/ROMSX_make_new_level.cpp +++ b/Source/ROMSX_make_new_level.cpp @@ -18,8 +18,8 @@ void ROMSX::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, const DistributionMapping& dm) { - cons_new[lev]->define(ba, dm, Cons::NumVars, cons_new[lev-1]->nGrowVect()); - cons_old[lev]->define(ba, dm, Cons::NumVars, cons_new[lev-1]->nGrowVect()); + cons_new[lev]->define(ba, dm, NCONS, cons_new[lev-1]->nGrowVect()); + cons_old[lev]->define(ba, dm, NCONS, cons_new[lev-1]->nGrowVect()); xvel_new[lev]->define(convert(ba, IntVect(1,0,0)), dm, 1, xvel_new[lev-1]->nGrowVect()); xvel_old[lev]->define(convert(ba, IntVect(1,0,0)), dm, 1, xvel_new[lev-1]->nGrowVect()); @@ -56,8 +56,8 @@ ROMSX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa int ngrow_vels = ComputeGhostCells(solverChoice.spatial_order)+2; #endif - MultiFab* tmp_cons_new; tmp_cons_new->define(ba, dm, Cons::NumVars, ngrow_state); - MultiFab* tmp_cons_old; tmp_cons_old->define(ba, dm, Cons::NumVars, ngrow_state); + MultiFab* tmp_cons_new; tmp_cons_new->define(ba, dm, NCONS, ngrow_state); + MultiFab* tmp_cons_old; tmp_cons_old->define(ba, dm, NCONS, ngrow_state); MultiFab* tmp_xvel_new; tmp_xvel_new->define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); MultiFab* tmp_xvel_old; tmp_xvel_old->define(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); @@ -74,7 +74,7 @@ ROMSX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa FillPatch(lev, time, tmp_yvel_new, yvel_new); FillPatch(lev, time, tmp_zvel_new, zvel_new); - MultiFab::Copy(*tmp_cons_old,*tmp_cons_new,0,0,NVAR,tmp_cons_new[lev].nGrowVect()); + MultiFab::Copy(*tmp_cons_old,*tmp_cons_new,0,0,NCONS,tmp_cons_new[lev].nGrowVect()); MultiFab::Copy(*tmp_xvel_old,*tmp_xvel_new,0,0, 1,tmp_xvel_new[lev].nGrowVect()); MultiFab::Copy(*tmp_yvel_old,*tmp_yvel_new,0,0, 1,tmp_yvel_new[lev].nGrowVect()); MultiFab::Copy(*tmp_zvel_old,*tmp_zvel_new,0,0, 1,tmp_zvel_new[lev].nGrowVect()); @@ -118,8 +118,8 @@ void ROMSX::MakeNewLevelFromScratch (int lev, Real /*time*/, const BoxArray& ba, int ngrow_vels = ComputeGhostCells(solverChoice.spatial_order)+2; #endif - cons_old[lev] = new MultiFab(ba, dm, Cons::NumVars, ngrow_state); - cons_new[lev] = new MultiFab(ba, dm, Cons::NumVars, ngrow_state); + cons_old[lev] = new MultiFab(ba, dm, NCONS, ngrow_state); + cons_new[lev] = new MultiFab(ba, dm, NCONS, ngrow_state); xvel_new[lev] = new MultiFab(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); xvel_old[lev] = new MultiFab(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels); @@ -153,8 +153,7 @@ void ROMSX::resize_stuff(int lev) vec_visc3d_r.resize(lev+1); vec_visc2_p.resize(lev+1); vec_visc2_r.resize(lev+1); - vec_diff2_salt.resize(lev+1); - vec_diff2_temp.resize(lev+1); + vec_diff2.resize(lev+1); vec_ru.resize(lev+1); vec_rv.resize(lev+1); vec_rufrc.resize(lev+1); @@ -239,8 +238,7 @@ void ROMSX::init_stuff(int lev, const BoxArray& ba, const DistributionMapping& d // check dimensionality vec_visc2_p[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0))); // harmonic viscosity at psi points -- difference to 3d? vec_visc2_r[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0))); // harmonic viscosity at rho points - vec_diff2_salt[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0))); // harmonic diffusivity salt - vec_diff2_temp[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0))); // harmonic diffusivity temperature + vec_diff2[lev].reset(new MultiFab(ba,dm,2,IntVect(NGROW,NGROW,0))); // harmonic diffusivity temperature/salt // maybe TODO: clean up component indexing in prestep? vec_ru[lev].reset(new MultiFab(convert(ba,IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,NGROW))); // RHS u (incl horizontal and vertical advection) diff --git a/Source/TimeIntegration/ROMSX_ComputeTimestep.cpp b/Source/TimeIntegration/ROMSX_ComputeTimestep.cpp index e7ab4ad9..c7f708a6 100644 --- a/Source/TimeIntegration/ROMSX_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ROMSX_ComputeTimestep.cpp @@ -64,7 +64,7 @@ ROMSX::estTimeStep(int level) const return new_lm_dt; }); - amrex::ParallelDescriptor::ReduceRealMax(estdt_lowM_inv); + ParallelDescriptor::ReduceRealMax(estdt_lowM_inv); if (estdt_lowM_inv > 0.0_rt) estdt_lowM = cfl / estdt_lowM_inv;; diff --git a/Source/TimeIntegration/ROMSX_advance_3d.cpp b/Source/TimeIntegration/ROMSX_advance_3d.cpp index c6ae98a5..cdb8f6df 100644 --- a/Source/TimeIntegration/ROMSX_advance_3d.cpp +++ b/Source/TimeIntegration/ROMSX_advance_3d.cpp @@ -25,12 +25,12 @@ ROMSX::advance_3d (int lev, std::unique_ptr& mf_Hvom, std::unique_ptr& mf_z_w, std::unique_ptr& mf_h, - const int ncomp, const int N, Real dt_lev) + const int N, Real dt_lev) { const int Mm = Geom(lev).Domain().size()[1]; - const int nrhs = ncomp-1; - const int nnew = ncomp-1; + const int nrhs = 0; + const int nnew = 0; const GpuArray dx = Geom(lev).CellSizeArray(); const GpuArray dxi = Geom(lev).InvCellSizeArray(); @@ -123,8 +123,7 @@ ROMSX::advance_3d (int lev, // // Update to u and v // - ParallelFor(amrex::makeSlab(tbxp2,2,0), - [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(amrex::makeSlab(tbxp2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { pm(i,j,0)=dxi[0]; pn(i,j,0)=dxi[1]; @@ -141,14 +140,14 @@ ROMSX::advance_3d (int lev, } else { cff=0.25*dt_lev*23.0/12.0; } - amrex::ParallelFor(xbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + + ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { u(i,j,k) += cff * (pm(i,j,0)+pm(i-1,j,0)) * (pn(i,j,0)+pn(i-1,j,0)) * ru(i,j,k,nrhs); u(i,j,k) *= 2.0 / (Hz(i-1,j,k) + Hz(i,j,k)); }); - amrex::ParallelFor(ybx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { v(i,j,k) += cff * (pm(i,j,0)+pm(i,j-1,0)) * (pn(i,j,0)+pn(i,j-1,0)) * rv(i,j,k,nrhs); v(i,j,k) *= 2.0 / (Hz(i,j-1,k) + Hz(i,j,k)); @@ -261,10 +260,10 @@ ROMSX::advance_3d (int lev, auto pm = fab_pm.array(); auto W= fab_W.array(); + // // Update to u and v // - ParallelFor(amrex::makeSlab(tbxp2,2,0), - [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(amrex::makeSlab(tbxp2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { pm(i,j,0)=dxi[0]; pn(i,j,0)=dxi[1]; @@ -310,8 +309,7 @@ ROMSX::advance_3d (int lev, } }); - ParallelFor(makeSlab(gbx1,2,N), - [=] AMREX_GPU_DEVICE (int i, int j, int) + ParallelFor(makeSlab(gbx1,2,N), [=] AMREX_GPU_DEVICE (int i, int j, int) { W(i,j,N) = 0.0; }); @@ -379,8 +377,7 @@ ROMSX::advance_3d (int lev, // // Update to u and v // - ParallelFor(makeSlab(tbxp2,2,0), - [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(makeSlab(tbxp2,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { pm(i,j,0)=dxi[0]; pn(i,j,0)=dxi[1]; diff --git a/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp b/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp index f3edf691..1ded00b2 100644 --- a/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp +++ b/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp @@ -13,16 +13,10 @@ void ROMSX::advance_3d_ml (int lev, Real dt_lev) FillPatch(lev, t_old[lev], zvel_old[lev], zvel_old); MultiFab mf_temp(*cons_new[lev], amrex::make_alias, Temp_comp, 1); -#ifdef ROMSX_USE_SALINITY MultiFab mf_salt(*cons_new[lev], amrex::make_alias, Salt_comp, 1); -#else - MultiFab mf_salt(*cons_new[lev], amrex::make_alias, Temp_comp, 1); -#endif auto N = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ - const int ncomp = 1; - advance_3d(lev, *xvel_new[lev], *yvel_new[lev], mf_temp, mf_salt, vec_t3[lev].get(), vec_s3[lev].get(), @@ -31,7 +25,7 @@ void ROMSX::advance_3d_ml (int lev, Real dt_lev) vec_DV_avg1[lev], vec_DV_avg2[lev], vec_ubar[lev], vec_vbar[lev], vec_Akv[lev], vec_Akt[lev], vec_Hz[lev], vec_Huon[lev], vec_Hvom[lev], - vec_z_w[lev], vec_hOfTheConfusingName[lev], ncomp, N, dt_lev); + vec_z_w[lev], vec_hOfTheConfusingName[lev], N, dt_lev); vec_ubar[lev]->FillBoundary(geom[lev].periodicity()); vec_vbar[lev]->FillBoundary(geom[lev].periodicity()); diff --git a/Source/TimeIntegration/ROMSX_coriolis.cpp b/Source/TimeIntegration/ROMSX_coriolis.cpp index 78a409ea..23fb4721 100644 --- a/Source/TimeIntegration/ROMSX_coriolis.cpp +++ b/Source/TimeIntegration/ROMSX_coriolis.cpp @@ -19,7 +19,7 @@ ROMSX::coriolis (const Box& xbx, const Box& ybx, //----------------------------------------------------------------------- // - amrex::ParallelFor(xbx, + ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real UFx_i = 0.5 * Hz(i ,j,k) * fomn(i ,j,0) * (vold(i ,j,k,nrhs)+vold(i ,j+1,k,nrhs)); @@ -27,7 +27,7 @@ ROMSX::coriolis (const Box& xbx, const Box& ybx, ru(i,j,k,nr) += 0.5*(UFx_i + UFx_im1); }); - amrex::ParallelFor(ybx, + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real VFe_j = 0.5 * Hz(i,j ,k) * fomn(i,j ,0) * (uold(i,j ,k,nrhs)+uold(i+1,j ,k,nrhs)); diff --git a/Source/TimeIntegration/ROMSX_prestep.cpp b/Source/TimeIntegration/ROMSX_prestep.cpp index a6fff0cc..977c40e0 100644 --- a/Source/TimeIntegration/ROMSX_prestep.cpp +++ b/Source/TimeIntegration/ROMSX_prestep.cpp @@ -12,10 +12,8 @@ using namespace amrex; * @param[inout] mf_v maybe just out * @param[inout] mf_ru * @param[inout] mf_rv - * @param[in ] mf_tempold - * @param[in ] mf_saltold - * @param[inout] mf_temp maybe just out - * @param[inout] mf_salt maybe just out + * @param[in ] S_old + * @param[inout] S_new * @param[out ] mf_W * @param[none ] mf_DC (temp) * @param[in ] mf_z_r @@ -40,8 +38,7 @@ ROMSX::prestep (int lev, MultiFab& mf_u, MultiFab& mf_v, std::unique_ptr& mf_ru, std::unique_ptr& mf_rv, - MultiFab& mf_tempold, MultiFab& mf_saltold, - MultiFab& mf_temp, MultiFab& mf_salt, + MultiFab& S_old, MultiFab& S_new, MultiFab& mf_W, MultiFab& mf_DC, std::unique_ptr& mf_z_r, std::unique_ptr& mf_z_w, @@ -59,17 +56,18 @@ ROMSX::prestep (int lev, const auto dx = Geom(lev).CellSizeArray(); const int Mm = Geom(lev).Domain().size()[1]; - const BoxArray& ba = mf_temp.boxArray(); - const DistributionMapping& dm = mf_temp.DistributionMap(); + const BoxArray& ba = S_old.boxArray(); + const DistributionMapping& dm = S_old.DistributionMap(); // Maybe not the best way to do this, but need to cache salt and temp since // they get rewritten by prestep_t MultiFab mf_saltcache(ba,dm,1,IntVect(NGROW,NGROW,0)); MultiFab mf_tempcache(ba,dm,1,IntVect(NGROW,NGROW,0)); - MultiFab::Copy(mf_saltcache,mf_salt,0,0,mf_salt.nComp(),IntVect(NGROW,NGROW,0)); //mf_salt.nGrowVect()); - MultiFab::Copy(mf_tempcache,mf_temp,0,0,mf_temp.nComp(),IntVect(NGROW,NGROW,0)); - for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + MultiFab::Copy(mf_saltcache,S_new,Salt_comp,0,1,IntVect(NGROW,NGROW,0)); + MultiFab::Copy(mf_tempcache,S_new,Temp_comp,0,1,IntVect(NGROW,NGROW,0)); + + for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& DC = mf_DC.array(mfi); Array4 const& Akv = (vec_Akv[lev])->array(mfi); @@ -84,10 +82,12 @@ ROMSX::prestep (int lev, Array4 const& vold = (mf_vold).array(mfi); Array4 const& u = (mf_u).array(mfi); Array4 const& v = (mf_v).array(mfi); - Array4 const& tempold = (mf_tempold).array(mfi); - Array4 const& saltold = (mf_saltold).array(mfi); - Array4 const& temp = (mf_temp).array(mfi); - Array4 const& salt = (mf_salt).array(mfi); + + Array4 const& tempold = S_old.array(mfi,Temp_comp); + Array4 const& saltold = S_old.array(mfi,Salt_comp); + Array4 const& temp = S_new.array(mfi,Temp_comp); + Array4 const& salt = S_new.array(mfi,Salt_comp); + Array4 const& tempstore = (vec_t3[lev])->array(mfi); Array4 const& saltstore = (vec_s3[lev])->array(mfi); Array4 const& ru = (mf_ru)->array(mfi); @@ -169,76 +169,29 @@ ROMSX::prestep (int lev, auto pnom_v=fab_pnom_v.array(); - amrex::ParallelFor(tbxp2D, + ParallelFor(tbxp2D, [=] AMREX_GPU_DEVICE (int i, int j, int ) { - //Note: are the comment definitons right? Don't seem to match metrics.f90 - pm(i,j,0) = dxi[0]; - pn(i,j,0) = dxi[1]; - om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j)) - on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j)) - om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j) - on_r(i,j,0)=1.0/dxi[1]; // 1/pn(i,j) - //todo: om_p on_p - om_p(i,j,0)=1.0/dxi[0]; // 4/(pm(i-1,j-1)+pm(i-1,j)+pm(i,j-1)+pm(i,j)) - on_p(i,j,0)=1.0/dxi[1]; // 4/(pn(i-1,j-1)+pn(i-1,j)+pn(i,j-1)+pn(i,j)) - on_v(i,j,0)=1.0/dxi[1]; // 2/(pn(i-1,j)+pn(i,j)) - om_u(i,j,0)=1.0/dxi[0]; // 2/(pm(i-1,j)+pm(i,j)) - pmon_u(i,j,0)=1.0; // (pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j)) - pnom_u(i,j,0)=1.0; // (pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j)) - pmon_v(i,j,0)=1.0; // (pm(i,j-1)+pm(i,j))/(pn(i,j-1)+pn(i,j)) - pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j)) + //Note: are the comment definitons right? Don't seem to match metrics.f90 + pm(i,j,0) = dxi[0]; + pn(i,j,0) = dxi[1]; }); - if (verbose > 2) { - amrex::PrintToFile("temp_preprestep").SetPrecision(18)< 1) { - Print() << "Akv box " << Box(Akv) << std::endl; - } prestep_t_3d(bx, gbx, tempold, temp, tempcache, ru, Hz, Akt, Huon, Hvom, pm, pn, W, DC, FC, tempstore, z_r, z_w, h, iic, ntfirst, nnew, nstp, nrhs, N, - lambda, dt_lev); + lambda, dt_lev); prestep_t_3d(bx, gbx, saltold, salt, saltcache, ru, Hz, Akt, Huon, Hvom, pm, pn, W, DC, FC, saltstore, z_r, z_w, h, iic, ntfirst, nnew, nstp, nrhs, N, - lambda, dt_lev); - - if (verbose > 2) { - amrex::PrintToFile("u").SetPrecision(18)< 2) { - amrex::PrintToFile("u_after_prestep").SetPrecision(18)< 1) { - amrex::AllPrint() << "Box(Huon) " << Box(Huon) << std::endl; - amrex::AllPrint() << "Box(Hvom) " << Box(Hvom) << std::endl; - } //Use FC and DC as intermediate arrays for FX and FE //First pass do centered 2d terms if (solverChoice.flat_bathymetry) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FX(i,j,k)=Box(tempold).contains(i-1,j,k) ? Huon(i,j,k)* 0.5*(tempold(i-1,j,k)+ tempold(i ,j,k)) : 1e34; }); - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FE(i,j,k)=Box(tempold).contains(i,j-1,k) ? Hvom(i,j,k)* @@ -162,13 +158,13 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, }); } else { - amrex::ParallelFor(utbxp1, + ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 FX(i,j,k)=tempold(i,j,k,nrhs)-tempold(i-1,j,k,nrhs); }); - amrex::ParallelFor(vtbxp1, + ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 @@ -178,14 +174,14 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, Real cffa=1.0/6.0; Real cffb=1.0/3.0; if (solverChoice.Hadv_scheme == AdvectionScheme::upstream3) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //Upstream3 curv(i,j,k)=-FX(i,j,k)+FX(i+1,j,k); }); //HACK to avoid using the wrong index of t (using upstream3) - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real max_Huon=max(Huon(i,j,k),0.0); //FArrayBox(Huon).max(); @@ -195,13 +191,13 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, }); } else if (solverChoice.Hadv_scheme == AdvectionScheme::centered4) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //Centered4 grad(i,j,k)=0.5*(FX(i,j,k)+FX(i+1,j,k)); }); - amrex::ParallelFor(ubx, + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FX(i,j,k)=Huon(i,j,k)*0.5*(tempold(i,j,k)+tempold(i-1,j,k)- @@ -212,13 +208,13 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, Error("Not a valid horizontal advection scheme"); } if (solverChoice.Hadv_scheme == AdvectionScheme::upstream3) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k); }); //HACK to avoid using the wrong index of t (using upstream3) - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real max_Hvom=max(Hvom(i,j,k),0.0); //FArrayBox(Huon).max(); @@ -228,12 +224,12 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, }); } else if (solverChoice.Hadv_scheme == AdvectionScheme::centered4) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { grad(i,j,k)=0.5*(FE(i,j,k)+FE(i,j+1,k)); }); - amrex::ParallelFor(vbx, + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FE(i,j,k)=Hvom(i,j,k)*0.5*(tempold(i,j,k)+tempold(i,j-1,k)- @@ -264,20 +260,8 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, cff1=0.5+GammaT; cff2=0.5-GammaT; } - if (verbose > 1) { - Print() << std::endl; - Print()<<(Box(tempstore))< 2) { - PrintToFile("prestep_before_t_tempstore").SetPrecision(18) << FArrayBox(tempstore) << std::endl; - PrintToFile("prestep_before_t_tempcache").SetPrecision(18) << FArrayBox(tempcache) << std::endl; - PrintToFile("prestep_before_t_tempold").SetPrecision(18) << FArrayBox(tempold) << std::endl; - } - amrex::ParallelFor(tbx, + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { tempstore(i,j,k)=Hz(i,j,k)*(cff1*tempold(i,j,k)+ @@ -292,63 +276,42 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, (FC(i+1,j)-FC(i,j)+ DC(i,j+1)-DC(i,j));*/ }); - if (verbose > 2) { - PrintToFile("prestep_t_tempstore").SetPrecision(18) << FArrayBox(tempstore) << std::endl; - PrintToFile("prestep_t_tempcache").SetPrecision(18) << FArrayBox(tempcache) << std::endl; - PrintToFile("prestep_t_tempold").SetPrecision(18) << FArrayBox(tempold) << std::endl; - } // // Time-step vertical advection of tracers (Tunits). Impose artificial // continuity equation. // - amrex::ParallelFor(tbx, + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //----------------------------------------------------------------------- // Add in vertical advection. //----------------------------------------------------------------------- - Real cff1=0.5; - Real cff2=7.0/12.0; - Real cff3=1.0/12.0; + Real c1=0.5; + Real c2=7.0/12.0; + Real c3=1.0/12.0; if (k>=1 && k<=N-2) { - FC(i,j,k)=( cff2*(tempold(i ,j,k ,nrhs)+ tempold(i,j,k+1,nrhs)) - -cff3*(tempold(i ,j,k-1,nrhs)+ tempold(i,j,k+2,nrhs)) )* + FC(i,j,k)=( c2*(tempold(i ,j,k ,nrhs)+ tempold(i,j,k+1,nrhs)) + -c3*(tempold(i ,j,k-1,nrhs)+ tempold(i,j,k+2,nrhs)) )* ( W(i,j,k)); } else // this needs to be split up so that the following can be concurrent { FC(i,j,N)=0.0; - FC(i,j,N-1)=( cff2*tempold(i ,j,N-1,nrhs)+ cff1*tempold(i,j,N ,nrhs) - -cff3*tempold(i ,j,N-2,nrhs) )* - ( W(i ,j,N-1)); + FC(i,j,N-1) = ( c2*tempold(i,j,N-1,nrhs)+ c1*tempold(i,j,N,nrhs)-c3*tempold(i,j,N-2,nrhs) ) + * W(i,j,N-1); - FC(i,j,0)=( cff2*tempold(i ,j,1,nrhs)+ cff1*tempold(i,j,0,nrhs) - -cff3*tempold(i ,j,2,nrhs) )* - ( W(i ,j,0)); - - // FC(i,0,-1)=0.0; + FC(i,j, 0) = ( c2*tempold(i,j, 1,nrhs)+ c1*tempold(i,j,0,nrhs)-c3*tempold(i,j,2,nrhs) ) + * W(i,j,0); } }); - // exit(0); -/* - Real cff; - - Real GammaT = 1.0/6.0; - if (iic==ntfirst) - { - cff=0.5*dt_lev; - } else { - cff=(1-GammaT)*dt_lev; - }*/ - - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { if(k-1>=0) { @@ -373,18 +336,14 @@ ROMSX::prestep_t_3d (const Box& tbx, const Box& gbx, DC(i,j+1)-DC(i,j));*/ }); - amrex::ParallelFor(tbx, + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real cff1=cff*pm(i,j,0)*pn(i,j,0); - Real cff4; - if(k-1>=0) { - cff4=FC(i,j,k)-FC(i,j,k-1); - } else { - cff4=FC(i,j,k); - } - tempstore(i,j,k)=DC(i,j,k)*(tempstore(i,j,k)-cff1*cff4); -// temp(i,j,k)=tempold(i,j,k); + Real c1 = cff*pm(i,j,0)*pn(i,j,0); + + Real c4 = (k>0) ? FC(i,j,k)-FC(i,j,k-1) : FC(i,j,k); + + tempstore(i,j,k) = DC(i,j,k)*(tempstore(i,j,k)-c1*c4); }); //----------------------------------------------------------------------- diff --git a/Source/TimeIntegration/ROMSX_prestep_uv_3d.cpp b/Source/TimeIntegration/ROMSX_prestep_uv_3d.cpp index a84071a1..e9e47a3a 100644 --- a/Source/TimeIntegration/ROMSX_prestep_uv_3d.cpp +++ b/Source/TimeIntegration/ROMSX_prestep_uv_3d.cpp @@ -12,9 +12,6 @@ ROMSX::prestep_uv_3d (const Box& tbx, const Box& gbx, Array4 u , Array4 v, Array4 ru, Array4 rv, Array4 Hz, Array4 Akv, - Array4 /*on_u*/, - Array4 /*om_v*/, - Array4 Huon, Array4 Hvom, Array4 pm, Array4 pn, Array4 /*W */, Array4 DC, Array4 FC , Array4 z_r, @@ -36,11 +33,6 @@ ROMSX::prestep_uv_3d (const Box& tbx, const Box& gbx, gbx2.grow(IntVect(NGROW,NGROW,0)); gbx1.grow(IntVect(NGROW-1,NGROW-1,0)); - if (verbose > 1) { - amrex::AllPrint() << "Box(Huon) " << Box(Huon) << std::endl; - amrex::AllPrint() << "Box(Hvom) " << Box(Hvom) << std::endl; - } - //Need to include pre_step3d.F terms update_vel_3d(tbxp1, gbx, 1, 0, u, uold, ru, Hz, Akv, DC, FC, diff --git a/Source/TimeIntegration/ROMSX_prsgrd.cpp b/Source/TimeIntegration/ROMSX_prsgrd.cpp index c4c42662..4b5e00e3 100644 --- a/Source/TimeIntegration/ROMSX_prsgrd.cpp +++ b/Source/TimeIntegration/ROMSX_prsgrd.cpp @@ -48,7 +48,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, auto dZx=fab_dZx.array(); // Derivatives in the z direction - amrex::ParallelFor(phi_bx, + ParallelFor(phi_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { if(k>=0&&k=0;k--) { @@ -78,7 +78,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, } }); - amrex::ParallelFor(phi_bxD, + ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) { Real cff1=1.0_rt/(z_r(i,j,N)-z_r(i,j,N-1)); @@ -100,7 +100,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, //This should be nodal // Derivatives in the x direction - amrex::ParallelFor(phi_ubx, + ParallelFor(phi_ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FC(i,j,k)=rho(i,j,k)-rho(i-1,j,k); @@ -109,7 +109,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, //This should be nodal aux and FC need wider boxes above //dZx and dRx may have index mismatch issues at k=2 and k=N - amrex::ParallelFor(phi_bxD, + ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) { for(int k=N;k>=0;k--) { @@ -131,7 +131,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, }); //This should be nodal aux and FC need wider boxes above - amrex::ParallelFor(utbxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(utbxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) { for(int k=N;k>=0;k--) { @@ -150,7 +150,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, }); //This should be nodal - amrex::ParallelFor(phi_vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(phi_vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FC(i,j,k)= rho(i,j,k)-rho(i,j-1,k); aux(i,j,k)= z_r(i,j,k)-z_r(i,j-1,k); @@ -158,7 +158,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, //This should be nodal aux and FC need wider boxes above //dZx and dRx may have index mismatch issues at k=2 and k=N - amrex::ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) { for(int k=N;k>=0;k--) { Real cff= 2.0*aux(i,j,k)*aux(i,j+1,k); @@ -179,7 +179,7 @@ ROMSX::prsgrd (const Box& phi_bx, const Box& phi_gbx, }); //This should be nodal aux and FC need wider boxes above - amrex::ParallelFor(vtbxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(vtbxD, [=] AMREX_GPU_DEVICE (int i, int j, int ) { for (int k=N;k>=0;k--) { diff --git a/Source/TimeIntegration/ROMSX_rho_eos.cpp b/Source/TimeIntegration/ROMSX_rho_eos.cpp index 8385ab97..f4022401 100644 --- a/Source/TimeIntegration/ROMSX_rho_eos.cpp +++ b/Source/TimeIntegration/ROMSX_rho_eos.cpp @@ -6,28 +6,25 @@ using namespace amrex; * rho_eos * * @param[in ] bx box for calculation - * @param[in ] temp temperature + * @param[in ] state state holds temp, salt * @param[out] rho density * @param[out] rhoA vertically-averaged density * @param[out] rhoS density perturbation * @param[in ] Hz * @param[in ] z_w * @param[in ] h - * @param[in ] nrhs * @param[in ] N */ void ROMSX::rho_eos (const Box& bx, - Array4 temp, - Array4 salt, + Array4 state, Array4 rho, Array4 rhoA, Array4 rhoS, Array4 Hz, Array4 z_w, Array4 h, - const int nrhs, const int N) { // @@ -48,11 +45,12 @@ ROMSX::rho_eos (const Box& bx, Real T0 = solverChoice.T0; Real Tcoef = solverChoice.Tcoef; Real Scoef = solverChoice.Scoef; - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - rho(i,j,k) = R0 - R0*Tcoef*(temp(i,j,k,nrhs)-T0); - rho(i,j,k) += R0*Scoef*(salt(i,j,k,nrhs)-S0) - 1000.0_rt; + rho(i,j,k) = R0 - R0*Tcoef*(state(i,j,k,Temp_comp)-T0) + + R0*Scoef*(state(i,j,k,Salt_comp)-S0) + - 1000.0_rt; }); // @@ -63,7 +61,7 @@ ROMSX::rho_eos (const Box& bx, // Real cff2 =1.0_rt/solverChoice.rho0; - amrex::ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { Real cff0 = rho(i,j,N)*Hz(i,j,N); rhoS(i,j,0) = 0.5_rt*cff0*Hz(i,j,N); diff --git a/Source/TimeIntegration/ROMSX_rhs_t_3d.cpp b/Source/TimeIntegration/ROMSX_rhs_t_3d.cpp index 372f3b31..1b522700 100644 --- a/Source/TimeIntegration/ROMSX_rhs_t_3d.cpp +++ b/Source/TimeIntegration/ROMSX_rhs_t_3d.cpp @@ -71,13 +71,13 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, fab_FX.template setVal(0.); fab_FE.template setVal(0.); - amrex::ParallelFor(bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { oHz(i,j,k) = 1.0/ Hz(i,j,k); }); - amrex::ParallelFor(utbxp1, + ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 @@ -91,7 +91,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, if (solverChoice.Hadv_scheme == AdvectionScheme::upstream3) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //Upstream3 @@ -100,7 +100,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, //HACK to avoid using the wrong index of t (using upstream3) Real max_Huon=FArrayBox(Huon).max(); Real min_Huon=FArrayBox(Huon).min(); - amrex::ParallelFor(ubx, + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FX(i,j,k)=Huon(i,j,k)*0.5*(tempstore(i,j,k)+tempstore(i-1,j,k))+ @@ -109,13 +109,13 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } else if (solverChoice.Hadv_scheme == AdvectionScheme::centered4) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //Centered4 grad(i,j,k)=0.5*(FX(i,j,k)+FX(i+1,j,k)); }); - amrex::ParallelFor(ubx, + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FX(i,j,k)=Huon(i,j,k)*0.5*(tempstore(i,j,k)+tempstore(i-1,j,k))+ @@ -129,14 +129,14 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } else { if (solverChoice.Hadv_scheme == AdvectionScheme::upstream3) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //Upstream3 curv(i,j,k)=-FX(i,j,k)+FX(i+1,j,k); }); //HACK to avoid using the wrong index of t (using upstream3) - amrex::ParallelFor(ubx, + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real max_Huon=max(Huon(i,j,k),0.0); //FArrayBox(Huon).max(); @@ -147,13 +147,13 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } else if (solverChoice.Hadv_scheme == AdvectionScheme::centered4) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //Centered4 grad(i,j,k)=0.5*(FX(i,j,k)+FX(i+1,j,k)); }); - amrex::ParallelFor(ubx, + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FX(i,j,k)=Huon(i,j,k)*0.5*(tempstore(i,j,k)+tempstore(i-1,j,k)- @@ -165,7 +165,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } } // flat bathymetry? - amrex::ParallelFor(vtbxp1, + ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //should be t index 3 @@ -176,7 +176,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, cffb=1.0/3.0; if (solverChoice.flat_bathymetry) { if (solverChoice.Hadv_scheme == AdvectionScheme::upstream3) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k); @@ -184,7 +184,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, //HACK to avoid using the wrong index of t (using upstream3) Real max_Hvom=FArrayBox(Hvom).max(); Real min_Hvom=FArrayBox(Hvom).min(); - amrex::ParallelFor(vbx, + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FE(i,j,k)=Hvom(i,j,k)*0.5*(tempstore(i,j,k)+tempstore(i,j-1,k))+ @@ -192,12 +192,12 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, }); } else if (solverChoice.Hadv_scheme == AdvectionScheme::centered4) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { grad(i,j,k)=0.5*(FE(i,j,k)+FE(i,j+1,k)); }); - amrex::ParallelFor(vbx, + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FE(i,j,k)=Hvom(i,j,k)*0.5*(tempstore(i,j,k)+tempstore(i,j-1,k))+ @@ -210,13 +210,13 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } else { if (solverChoice.Hadv_scheme == AdvectionScheme::upstream3) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k); }); //HACK to avoid using the wrong index of t (using upstream3) - amrex::ParallelFor(vbx, + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real max_Hvom=max(Hvom(i,j,k),0.0); //FArrayBox(Huon).max(); @@ -226,12 +226,12 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, }); } else if (solverChoice.Hadv_scheme == AdvectionScheme::centered4) { - amrex::ParallelFor(tbxp1, + ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k) { grad(i,j,k)=0.5*(FE(i,j,k)+FE(i,j+1,k)); }); - amrex::ParallelFor(vbx, + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FE(i,j,k)=Hvom(i,j,k)*0.5*(tempstore(i,j,k)+tempstore(i,j-1,k)- @@ -243,7 +243,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } } - amrex::ParallelFor(bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // @@ -265,7 +265,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, // Fourth-order, central differences vertical advective flux // (Tunits m3/s). // - amrex::ParallelFor(bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { //----------------------------------------------------------------------- @@ -293,7 +293,7 @@ ROMSX::rhs_t_3d (const Box& bx, const Box& gbx, } }); - amrex::ParallelFor(bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real cff1=dt_lev*pm(i,j,0)*pn(i,j,0); diff --git a/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp b/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp index 954f20e7..c6f816d2 100644 --- a/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp +++ b/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp @@ -94,7 +94,7 @@ ROMSX::rhs_uv_2d (const Box& xbx, const Box& ybx, // // Add in horizontal advection. // - amrex::ParallelFor(makeSlab(xbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int) + ParallelFor(makeSlab(xbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int) { Real cff1=UFx(i,j ,0)-UFx(i-1,j,0); Real cff2=UFe(i,j+1,0)-UFe(i ,j,0); @@ -117,7 +117,7 @@ ROMSX::rhs_uv_2d (const Box& xbx, const Box& ybx, // and requires DUon on [ 0:nx-1,-2:ny+1] (0,0,0) y-faces // Grow ybx by one in high x-direction - amrex::ParallelFor(growHi(ybx,0,1), + ParallelFor(growHi(ybx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int ) { Real cff=1.0/6.0; @@ -138,7 +138,7 @@ ROMSX::rhs_uv_2d (const Box& xbx, const Box& ybx, VFx(i,j,0)=0.25*(cff1-vxx_avg*cff)* (cff2-cff*(Huee_j+ Huee_jm1)); }); - amrex::ParallelFor(growLo(ybx,1,1), + ParallelFor(growLo(ybx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int) { Real vee_j = vbar(i,j-1,0,krhs)-2.0*vbar(i,j ,0,krhs)+vbar(i,j+1,0,krhs); @@ -155,7 +155,7 @@ ROMSX::rhs_uv_2d (const Box& xbx, const Box& ybx, cff * (Hvee_j+ Hvee_jp1)); }); - amrex::ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int) + ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int) { Real cff1=VFx(i+1,j,0)-VFx(i ,j,0); Real cff2=VFe(i ,j,0)-VFe(i,j-1,0); diff --git a/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp b/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp index a3ec8f55..2f9a6305 100644 --- a/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp +++ b/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp @@ -11,13 +11,12 @@ using namespace amrex; * @param[in ] v velocity in y-direction * @param[out] Hv H-weighted velocity in y-direction * @param[in ] om_v weighting for v - * @param[in ] nrhs component of u and v */ void ROMSX::set_massflux_3d (Array4 u, Array4 Hu, Array4 on_u, Array4 v, Array4 Hv, Array4 om_v, - Array4 Hz, const int nrhs) + Array4 Hz) { // //----------------------------------------------------------------------- @@ -26,11 +25,11 @@ ROMSX::set_massflux_3d (Array4 u, Array4 Hu, Array4 on_u, // ParallelFor(Box(Hu), [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Hu(i,j,k)=0.5*(Hz(i,j,k)+Hz(i-1,j,k))*u(i,j,k,nrhs)* on_u(i,j,0); + Hu(i,j,k)=0.5*(Hz(i,j,k)+Hz(i-1,j,k))*u(i,j,k)* on_u(i,j,0); }); ParallelFor(Box(Hv), [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Hv(i,j,k)=0.5*(Hz(i,j,k)+Hz(i,j-1,k))*v(i,j,k,nrhs)* om_v(i,j,0); + Hv(i,j,k)=0.5*(Hz(i,j,k)+Hz(i,j-1,k))*v(i,j,k)* om_v(i,j,0); }); } diff --git a/Source/TimeIntegration/ROMSX_setup_step.cpp b/Source/TimeIntegration/ROMSX_setup_step.cpp index 41733938..ca713709 100644 --- a/Source/TimeIntegration/ROMSX_setup_step.cpp +++ b/Source/TimeIntegration/ROMSX_setup_step.cpp @@ -32,10 +32,8 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) const BoxArray& ba = S_old.boxArray(); const DistributionMapping& dm = S_old.DistributionMap(); - const int ncomp = 1; - const int nrhs = ncomp-1; - const int nnew = ncomp-1; - const int nstp = ncomp-1; + const int nrhs = 0; + const int nstp = 0; // Place-holder for source array -- for now just set to 0 MultiFab source(ba,dm,nvars,1); @@ -73,25 +71,11 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) std::unique_ptr& mf_rdrag = vec_rdrag[lev]; std::unique_ptr& mf_bustr = vec_bustr[lev]; std::unique_ptr& mf_bvstr = vec_bvstr[lev]; - MultiFab mf_temp(S_new, amrex::make_alias, Temp_comp, 1); -#ifdef ROMSX_USE_SALINITY - MultiFab mf_salt(S_new, amrex::make_alias, Salt_comp, 1); -#else - MultiFab mf_salt(S_new, amrex::make_alias, Temp_comp, 1); -#endif - MultiFab mf_tempold(S_old, amrex::make_alias, Temp_comp, 1); -#ifdef ROMSX_USE_SALINITY - MultiFab mf_saltold(S_old, amrex::make_alias, Salt_comp, 1); -#else - MultiFab mf_saltold(S_old, amrex::make_alias, Temp_comp, 1); -#endif MultiFab mf_rw(ba,dm,1,IntVect(NGROW,NGROW,0)); MultiFab mf_W(ba,dm,1,IntVect(NGROW+1,NGROW+1,0)); std::unique_ptr& mf_visc2_p = vec_visc2_p[lev]; std::unique_ptr& mf_visc2_r = vec_visc2_r[lev]; - std::unique_ptr& mf_diff2_temp = vec_diff2_temp[lev]; - std::unique_ptr& mf_diff2_salt = vec_diff2_salt[lev]; // We need to set these because otherwise in the first call to romsx_advance we may // read uninitialized data on ghost values in setting the bc's on the velocities @@ -116,10 +100,9 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) mf_vold.FillBoundary(geom[lev].periodicity()); mf_w.FillBoundary(geom[lev].periodicity()); mf_W.FillBoundary(geom[lev].periodicity()); - mf_tempold.FillBoundary(geom[lev].periodicity()); - mf_temp.FillBoundary(geom[lev].periodicity()); - mf_saltold.FillBoundary(geom[lev].periodicity()); - mf_salt.FillBoundary(geom[lev].periodicity()); + + S_old.FillBoundary(geom[lev].periodicity()); + S_new.FillBoundary(geom[lev].periodicity()); mf_rw.setVal(0.0); mf_W.setVal(0.0); @@ -146,7 +129,7 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) const int Mm = Geom(lev).Domain().size()[1]; //MFIter::allowMultipleMFIters(true); - for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& h = (vec_hOfTheConfusingName[lev])->const_array(mfi); Array4 const& Hz = (vec_Hz[lev])->const_array(mfi); @@ -158,8 +141,6 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) Array4 const& rho = (mf_rho).array(mfi); Array4 const& rhoA = (mf_rhoA)->array(mfi); Array4 const& rhoS = (mf_rhoS)->array(mfi); - Array4 const& tempold = (mf_tempold).const_array(mfi); - Array4 const& saltold = (mf_saltold).const_array(mfi); Array4 const& rdrag = (mf_rdrag)->array(mfi); Array4 const& bustr = (mf_bustr)->array(mfi); Array4 const& bvstr = (mf_bvstr)->array(mfi); @@ -222,16 +203,17 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) auto pmon_v=fab_pmon_v.array(); auto pnom_v=fab_pnom_v.array(); - amrex::ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int ) { on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j)) }); - amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int ) + + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int ) { om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j)) }); - amrex::ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (int i, int j, int ) { //Note: are the comment definitons right? Don't seem to match metrics.f90 om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j) @@ -247,33 +229,32 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j)) }); - amrex::ParallelFor(gbx2, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(gbx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Huon(i,j,k,0)=0.0; Hvom(i,j,k,0)=0.0; }); // Set bottom stress as defined in set_vbx.F - amrex::ParallelFor(gbx1D, - [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int ) { - bustr(i,j,0) = 0.5 * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0,nrhs)); - bvstr(i,j,0) = 0.5 * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0,nrhs)); + bustr(i,j,0) = 0.5 * (rdrag(i-1,j,0)+rdrag(i,j,0))*(uold(i,j,0)); + bvstr(i,j,0) = 0.5 * (rdrag(i,j-1,0)+rdrag(i,j,0))*(vold(i,j,0)); }); // Updates Huon and Hvom - set_massflux_3d(uold,Huon,on_u,vold,Hvom,om_v,Hz,nnew); + set_massflux_3d(uold,Huon,on_u,vold,Hvom,om_v,Hz); - rho_eos(gbx2,tempold,saltold,rho,rhoA,rhoS,Hz,z_w,h,nrhs,N); + Array4 const& state_old = S_old.const_array(mfi); + rho_eos(gbx2,state_old,rho,rhoA,rhoS,Hz,z_w,h,N); } if(solverChoice.use_prestep) { + const int nnew = 0; prestep(lev, mf_uold, mf_vold, - mf_u, mf_v, - mf_ru, mf_rv, mf_tempold, mf_saltold, - mf_temp, mf_salt, mf_W, + mf_u, mf_v, mf_ru, mf_rv, + S_old, S_new, mf_W, mf_DC, mf_z_r, mf_z_w, mf_h, mf_sustr, mf_svstr, mf_bustr, mf_bvstr, iic, ntfirst, nnew, nstp, nrhs, N, dt_lev); } @@ -281,7 +262,7 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) mf_W.FillBoundary(geom[lev].periodicity()); - for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + for ( MFIter mfi(S_old, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& Hz = (vec_Hz[lev])->array(mfi); Array4 const& Huon = (vec_Huon[lev])->array(mfi); @@ -293,8 +274,6 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) Array4 const& u = (mf_u).array(mfi); Array4 const& v = (mf_v).array(mfi); Array4 const& rho = (mf_rho).array(mfi); - Array4 const& temp = (mf_temp).array(mfi); - Array4 const& salt = (mf_salt).array(mfi); Array4 const& ru = (mf_ru)->array(mfi); Array4 const& rv = (mf_rv)->array(mfi); Array4 const& rufrc = (mf_rufrc)->array(mfi); @@ -306,8 +285,6 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) Array4 const& bvstr = (mf_bvstr)->array(mfi); Array4 const& visc2_p = (mf_visc2_p)->array(mfi); Array4 const& visc2_r = (mf_visc2_r)->array(mfi); - Array4 const& diff2_salt = (mf_diff2_salt)->array(mfi); - Array4 const& diff2_temp = (mf_diff2_temp)->array(mfi); Array4 const& zeta = (vec_zeta[lev])->array(mfi); Array4 const& Zt_avg1 = (vec_Zt_avg1[lev])->array(mfi); @@ -386,11 +363,11 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) Real coriolis_f0 = solverChoice.coriolis_f0; Real coriolis_beta = solverChoice.coriolis_beta; + const auto dx = Geom(lev).CellSizeArray(); + ParallelFor(tbxp2D, [=] AMREX_GPU_DEVICE (int i, int j, int ) { - const auto dx = geom[lev].CellSize(); - pm(i,j,0) = dxi[0]; pn(i,j,0) = dxi[1]; Real Esize=geom[lev].ProbHi()[1] - geom[lev].ProbLo()[1]; @@ -398,17 +375,17 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) Real f=coriolis_f0 + coriolis_beta*(y-.5*Esize); fomn(i,j,0)=f*(1.0/(pm(i,j,0)*pn(i,j,0))); }); - amrex::ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int ) { on_u(i,j,0)=1.0/dxi[1]; // 2/(pm(i,j-1)+pm(i,j)) }); - amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int ) + + ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int ) { om_v(i,j,0)=1.0/dxi[0]; // 2/(pm(i,j-1)+pm(i,j)) }); - amrex::ParallelFor(tbxp2D, - [=] AMREX_GPU_DEVICE (int i, int j, int ) + ParallelFor(tbxp2D, [=] AMREX_GPU_DEVICE (int i, int j, int ) { //Note: are the comment definitons right? Don't seem to match metrics.f90 om_r(i,j,0)=1.0/dxi[0]; // 1/pm(i,j) @@ -424,16 +401,18 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) pnom_v(i,j,0)=1.0; // (pn(i,j-1)+pn(i,j))/(pm(i,j-1)+pm(i,j)) }); - amrex::ParallelFor(gbx2, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(gbx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) { FC(i,j,k)=0.0; }); prsgrd(tbxp1,gbx1,utbx,vtbx,ru,rv,on_u,om_v,rho,FC,Hz,z_r,z_w,nrhs,N); - t3dmix(bx, temp, diff2_temp, Hz, pm, pn, pmon_u, pnom_v, nrhs, nnew, dt_lev); - t3dmix(bx, salt, diff2_salt, Hz, pm, pn, pmon_u, pnom_v, nrhs, nnew, dt_lev); + // Apply mixing to temperature and, if use_salt, salt + int ncomp = solverChoice.use_salt ? 2 : 1; + Array4 const& s_arr = S_old.array(mfi); + Array4 const& diff2_arr = vec_diff2[lev]->array(mfi); + t3dmix(bx, s_arr, diff2_arr, Hz, pm, pn, pmon_u, pnom_v, dt_lev, ncomp); if (solverChoice.use_coriolis) { //----------------------------------------------------------------------- @@ -457,12 +436,12 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) on_u, om_v, om_u, on_v, W, FC, nrhs, N); if(solverChoice.use_uv3dmix) { + const int nnew = 0; uv3dmix(xbx, ybx, u, v, uold, vold, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev); } // Set first two components of zeta to time-averaged values before barotropic update - amrex::ParallelFor(gbx2D, - [=] AMREX_GPU_DEVICE (int i, int j, int) + ParallelFor(gbx2D, [=] AMREX_GPU_DEVICE (int i, int j, int) { zeta(i,j,0,0) = Zt_avg1(i,j,0); zeta(i,j,0,1) = Zt_avg1(i,j,0); @@ -472,11 +451,8 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) // Update Akv with new depth. NOTE: this happens before set_zeta in ROMS set_vmix(lev); - mf_temp.FillBoundary(geom[lev].periodicity()); - mf_salt.FillBoundary(geom[lev].periodicity()); - - mf_tempold.FillBoundary(geom[lev].periodicity()); - mf_saltold.FillBoundary(geom[lev].periodicity()); + S_old.FillBoundary(geom[lev].periodicity()); + S_new.FillBoundary(geom[lev].periodicity()); vec_t3[lev]->FillBoundary(geom[lev].periodicity()); vec_s3[lev]->FillBoundary(geom[lev].periodicity()); diff --git a/Source/TimeIntegration/ROMSX_t3dmix.cpp b/Source/TimeIntegration/ROMSX_t3dmix.cpp index f1c688a1..226017ab 100644 --- a/Source/TimeIntegration/ROMSX_t3dmix.cpp +++ b/Source/TimeIntegration/ROMSX_t3dmix.cpp @@ -4,12 +4,11 @@ using namespace amrex; void ROMSX::t3dmix (const Box& bx, - Array4 t, + Array4 state, Array4 diff2, Array4 Hz, Array4 pm, Array4 pn, Array4 pmon_u, Array4 pnom_v, - int nrhs, int nnew, - const amrex::Real dt_lev) + const Real dt_lev, const int ncomp) { //----------------------------------------------------------------------- // Add in harmonic diffusivity s terms. @@ -18,40 +17,32 @@ ROMSX::t3dmix (const Box& bx, Box xbx(bx); xbx.surroundingNodes(0); Box ybx(bx); ybx.surroundingNodes(1); - FArrayBox fab_FX(xbx,1,amrex::The_Async_Arena()); - FArrayBox fab_FE(ybx,1,amrex::The_Async_Arena()); - - fab_FX.template setVal(0.); - fab_FE.template setVal(0.); + FArrayBox fab_FX(xbx,ncomp,The_Async_Arena()); + FArrayBox fab_FE(ybx,ncomp,The_Async_Arena()); auto FX=fab_FX.array(); auto FE=fab_FE.array(); - ParallelFor(xbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(xbx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - const amrex::Real cff = 0.25 * (diff2(i,j,0) + diff2(i-1,j,0)) * pmon_u(i,j,0); - FX(i,j,k) = cff * (Hz(i,j,k)+Hz(i+1,j,k))*(t(i,j,k,nrhs)-t(i-1,j,k,nrhs)); + const Real cff = 0.25 * (diff2(i,j,n) + diff2(i-1,j,n)) * pmon_u(i,j,0); + FX(i,j,k,n) = cff * (Hz(i,j,k) + Hz(i+1,j,k)) * (state(i,j,k,n)-state(i-1,j,k,n)); }); - ParallelFor(ybx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(ybx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - const amrex::Real cff=0.25*(diff2(i,j,0)+diff2(i,j-1,0)) * pnom_v(i,j,0); - FE(i,j,k) = cff * (Hz(i,j,k) + Hz(i,j-1,k)) * (t(i,j,k,nrhs) - t(i,j-1,k,nrhs)); + const Real cff = 0.25*(diff2(i,j,n)+diff2(i,j-1,n)) * pnom_v(i,j,0); + FE(i,j,k,n) = cff * (Hz(i,j,k) + Hz(i,j-1,k)) * (state(i,j,k,n) - state(i,j-1,k,n)); }); /* Time-step harmonic, S-surfaces diffusion term. */ - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { - const amrex::Real cff=dt_lev*pm(i,j,0)*pn(i,j,0); - const amrex::Real cff1=cff*(FX(i+1,j ,k)-FX(i,j,k)); - const amrex::Real cff2=cff*(FE(i ,j+1,k)-FE(i,j,k)); - const amrex::Real cff3=cff1+cff2; + const Real cff = dt_lev*pm(i,j,0)*pn(i,j,0); - t(i,j,k,nnew) += cff3; + state(i,j,k,n) += cff * ( (FX(i+1,j ,k,n)-FX(i,j,k,n)) + +(FE(i ,j+1,k,n)-FE(i,j,k,n)) ); }); } diff --git a/Source/TimeIntegration/ROMSX_update_vel_3d.cpp b/Source/TimeIntegration/ROMSX_update_vel_3d.cpp index af3363e1..dac308ff 100644 --- a/Source/TimeIntegration/ROMSX_update_vel_3d.cpp +++ b/Source/TimeIntegration/ROMSX_update_vel_3d.cpp @@ -49,7 +49,7 @@ ROMSX::update_vel_3d (const Box& vel_bx, const Box& gbx, Print() << "vel old " << Box(vel_old) << std::endl; Print() << "Akv " << Box(Akv) << std::endl; } - amrex::ParallelFor(vel_bx, + ParallelFor(vel_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { if(k+1<=N&&k>=0) @@ -72,7 +72,7 @@ ROMSX::update_vel_3d (const Box& vel_bx, const Box& gbx, * (pn(i,j,0)+pn(i-ioff,j-joff,0)); }); - amrex::ParallelFor(gbxvel, + ParallelFor(gbxvel, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real cff3=dt_lev*(1.0-lambda); diff --git a/Source/TimeIntegration/ROMSX_uv3dmix.cpp b/Source/TimeIntegration/ROMSX_uv3dmix.cpp index ad9d479b..481bc133 100644 --- a/Source/TimeIntegration/ROMSX_uv3dmix.cpp +++ b/Source/TimeIntegration/ROMSX_uv3dmix.cpp @@ -39,7 +39,7 @@ ROMSX::uv3dmix (const Box& xbx, const Box& ybx, // to do so requires UFe on [ 0:nx , 0:ny ] (1,1,0) xy-nodes // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces // and requires vold on [ -1:nx , 0:ny ] (0,1,0) y-faces - amrex::ParallelFor(growLo(xbx,0,1), + ParallelFor(growLo(xbx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { const amrex::Real cff = 0.5*Hz(i,j,k) * (pm(i,j,0) / pn(i,j,0) * @@ -50,7 +50,7 @@ ROMSX::uv3dmix (const Box& xbx, const Box& ybx, (pm(i,j-1,0)+pm(i,j ,0))*vold(i,j ,k,nrhs))); UFx(i,j,k) = on_r(i,j,0)*on_r(i,j,0)*visc2_r(i,j,0)*cff; }); - amrex::ParallelFor(growHi(xbx,1,1), + ParallelFor(growHi(xbx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { const amrex::Real cff = 0.125 * (Hz(i-1,j ,k)+Hz(i,j ,k)+ @@ -64,7 +64,7 @@ ROMSX::uv3dmix (const Box& xbx, const Box& ybx, UFe(i,j,k) = om_p(i,j,0)*om_p(i,j,0)*visc2_p(i,j,0)*cff; }); - amrex::ParallelFor(xbx, + ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { const amrex::Real cff=dt_lev*0.25*(pm(i-1,j,0)+pm(i,j,0))*(pn(i-1,j,0)+pn(i,j,0)); @@ -85,7 +85,7 @@ ROMSX::uv3dmix (const Box& xbx, const Box& ybx, // to do so requires VFx on [ 0:nx , 0:ny ] (1,1,0) xy-nodes // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces // and requires vold on [-1:nx , 0:ny ] (0,1,0) y-faces - amrex::ParallelFor(growLo(ybx,1,1), + ParallelFor(growLo(ybx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { // cff depends on k, but UFx and VFe will only be affected by the last cell? @@ -97,7 +97,7 @@ ROMSX::uv3dmix (const Box& xbx, const Box& ybx, (pm(i,j-1,0)+pm(i,j ,0))*vold(i,j ,k,nrhs))); VFe(i,j,k) = om_r(i,j,0)*om_r(i,j,0)*visc2_r(i,j,0)*cff; }); - amrex::ParallelFor(growHi(ybx,0,1), + ParallelFor(growHi(ybx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k) { const amrex::Real cff = 0.125 * (Hz(i-1,j ,k)+Hz(i,j ,k)+ @@ -111,7 +111,7 @@ ROMSX::uv3dmix (const Box& xbx, const Box& ybx, VFx(i,j,k) = on_p(i,j,0)*on_p(i,j,0)*visc2_p(i,j,0)*cff; }); - amrex::ParallelFor(ybx, + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { const amrex::Real cff=dt_lev*0.25*(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0)); diff --git a/Tests/CTestList.cmake b/Tests/CTestList.cmake index 07d55ab1..60370460 100644 --- a/Tests/CTestList.cmake +++ b/Tests/CTestList.cmake @@ -94,8 +94,10 @@ endfunction(add_test_u) #============================================================================= if(WIN32) add_test_r(Upwelling "Upwelling/*/upwelling.exe" "plt00010") + add_test_r(Seamount "Seamount/*/seamount.exe" "plt00010") else() add_test_r(Upwelling "Upwelling/upwelling" "plt00010") + add_test_r(Seamount "Seamount/seamount" "plt00010") endif() #============================================================================= # Performance tests diff --git a/Tests/ROMSX_Gold_Files/Seamount/Header b/Tests/ROMSX_Gold_Files/Seamount/Header new file mode 100644 index 00000000..a76220e7 --- /dev/null +++ b/Tests/ROMSX_Gold_Files/Seamount/Header @@ -0,0 +1,40 @@ +HyperCLaw-V1.1 +6 +temp +omega +salt +x_velocity +y_velocity +z_velocity +3 +600 +0 +0 0 -5000 +320000 320000 0 + +((0,0,0) (48,47,12) (0,0,0)) +10 +6530.6122448979595 6666.666666666667 384.61538461538464 +0 +0 +0 6 600 +10 +0 111020.40816326531 +0 160000 +-5000 0 +0 111020.40816326531 +160000 320000 +-5000 0 +111020.40816326531 215510.20408163266 +0 160000 +-5000 0 +111020.40816326531 215510.20408163266 +160000 320000 +-5000 0 +215510.20408163266 320000 +0 160000 +-5000 0 +215510.20408163266 320000 +160000 320000 +-5000 0 +Level_0/Cell diff --git a/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00000 b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00000 new file mode 100644 index 00000000..6b0793ed Binary files /dev/null and b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00000 differ diff --git a/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00001 b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00001 new file mode 100644 index 00000000..8ce805c1 Binary files /dev/null and b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00001 differ diff --git a/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00002 b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00002 new file mode 100644 index 00000000..eef623cf Binary files /dev/null and b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00002 differ diff --git a/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00003 b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00003 new file mode 100644 index 00000000..52e68ad1 Binary files /dev/null and b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_D_00003 differ diff --git a/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_H b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_H new file mode 100644 index 00000000..998afff2 --- /dev/null +++ b/Tests/ROMSX_Gold_Files/Seamount/Level_0/Cell_H @@ -0,0 +1,36 @@ +1 +1 +6 +0 +(6 0 +((0,0,0) (16,23,12) (0,0,0)) +((0,24,0) (16,47,12) (0,0,0)) +((17,0,0) (32,23,12) (0,0,0)) +((17,24,0) (32,47,12) (0,0,0)) +((33,0,0) (48,23,12) (0,0,0)) +((33,24,0) (48,47,12) (0,0,0)) +) +6 +FabOnDisk: Cell_D_00002 0 +FabOnDisk: Cell_D_00000 0 +FabOnDisk: Cell_D_00003 0 +FabOnDisk: Cell_D_00000 254682 +FabOnDisk: Cell_D_00001 0 +FabOnDisk: Cell_D_00001 239706 + +6,6 +1.0091152195473329e+01,0.0000000000000000e+00,3.1999999999999950e+01,-8.0596562274061335e-06,-4.7571516124719847e-06,0.0000000000000000e+00, +1.0091152195473322e+01,0.0000000000000000e+00,3.1999999999999947e+01,-8.2337553238976310e-06,-7.3240047089317704e-06,0.0000000000000000e+00, +1.0091152209737174e+01,0.0000000000000000e+00,3.1999999999999943e+01,-1.7154954710763102e-04,-1.6816946567474648e-04,0.0000000000000000e+00, +1.0091152209737173e+01,0.0000000000000000e+00,3.1999999999999940e+01,-1.7555807622884347e-04,-1.0041814688398989e-04,0.0000000000000000e+00, +1.0091152195473322e+01,0.0000000000000000e+00,3.1999999999999947e+01,-1.4407271866654339e-05,-3.0824026304662890e-06,0.0000000000000000e+00, +1.0091152195473329e+01,0.0000000000000000e+00,3.1999999999999947e+01,-1.4375289106788807e-05,-6.1547967524413109e-06,0.0000000000000000e+00, + +6,6 +1.7436210557246753e+01,0.0000000000000000e+00,3.2000000000000050e+01,1.4375289102152935e-05,7.8067921994399393e-06,0.0000000000000000e+00, +1.7436210557231885e+01,0.0000000000000000e+00,3.2000000000000050e+01,1.4407271871353746e-05,4.4587673291830429e-06,0.0000000000000000e+00, +1.7464062300273408e+01,0.0000000000000000e+00,3.2000000000000057e+01,1.7555807622686217e-04,1.0041814688891416e-04,0.0000000000000000e+00, +1.7464062300273408e+01,0.0000000000000000e+00,3.2000000000000050e+01,1.7154954711114776e-04,1.6816946567141416e-04,0.0000000000000000e+00, +1.7434501861032235e+01,0.0000000000000000e+00,3.2000000000000050e+01,8.2337553251508738e-06,5.6784458386967229e-06,0.0000000000000000e+00, +1.7434501861035898e+01,0.0000000000000000e+00,3.2000000000000050e+01,8.0596562265545909e-06,3.3850743510283830e-06,0.0000000000000000e+00, + diff --git a/Tests/ROMSX_Gold_Files/Seamount/job_info b/Tests/ROMSX_Gold_Files/Seamount/job_info new file mode 100644 index 00000000..c3c70832 --- /dev/null +++ b/Tests/ROMSX_Gold_Files/Seamount/job_info @@ -0,0 +1,161 @@ +=============================================================================== + ROMSX Job Information +=============================================================================== +inputs file: inputs + +number of MPI processes: 4 + +CPU time used since start of simulation (CPU-hours): 0.0158847 + +=============================================================================== + Plotfile Information +=============================================================================== +output data / time: Tue Dec 12 10:35:43 2023 +output dir: /home/almgren/GitCode/ROMSX/Exec/Seamount + + +=============================================================================== + Build Information +=============================================================================== +build date: 2023-12-12 10:31:21.485092 +build machine: Linux manda 5.15.0-76-generic #83~20.04.1-Ubuntu SMP Wed Jun 21 20:23:31 UTC 2023 x86_64 x86_64 x86_64 GNU/Linux +build dir: /home/almgren/GitCode/ROMSX/Exec/Seamount +AMReX dir: /home/almgren/GitCode/amrex + +COMP: gnu +COMP version: 9.4.0 + + +ROMSX git hash: fe8d901-dirty +AMReX git hash: 23.12-8-g43d71da32-dirty + + +=============================================================================== + Grid Information +=============================================================================== + level: 0 + number of boxes = 6 + maximum zones = 49 48 13 + + Boundary conditions + -x: Periodic + +x: Periodic + -y: Periodic + +y: Periodic + -z: SlipWall + +z: SlipWall + + +=============================================================================== + Inputs File Parameters +=============================================================================== +max_step = 10 +stop_time = 30000.0 +amrex.fpe_trap_invalid = 1 +fabarray.mfiter_tile_size = 1024 1024 1024 +geometry.prob_lo = 0. 0. -5000. +geometry.prob_hi = 320000. 320000. 0. +amr.n_cell = 49 48 13 +geometry.is_periodic = 1 1 0 +zlo.type = SlipWall +zhi.type = SlipWall +romsx.fixed_dt = 60.0 +romsx.fixed_ndtfast_ratio = 20 +romsx.flat_bathymetry = 0 +romsx.sum_interval = 1 +romsx.v = 0 +amr.v = 1 +amr.max_level = 0 +romsx.check_file = chk +romsx.check_int = -57600 +romsx.plot_file_1 = plt +romsx.plot_int_1 = 1 +romsx.plot_vars_1 = omega salt temp x_velocity y_velocity z_velocity +romsx.plotfile_type = amrex +romsx.use_coriolis = true +romsx.horizontal_advection_scheme = upstream3 +romsx.spatial_order = 2 +prob.R0 = 1027.0 +prob.S0 = 32.0 +prob.T0 = 10.0 +prob.rho_0 = 1.0 +prob.T_0 = 1.0 +prob.A_0 = 1.0 +prob.u_0 = 0.0 +prob.v_0 = 0.0 +prob.rad_0 = 0.25 +prob.z0 = 0.1 +prob.zRef = 80.0e-3 +prob.uRef = 8.0e-3 +prob.prob_type = 12 +amr.refine_grid_layout_x = 1 +amr.refine_grid_layout_y = 1 +amr.refine_grid_layout_z = 0 +amr.n_proper = 2 +amr.max_grid_size = 2048 +amr.blocking_factor = 1 +amr.n_error_buf = 0 +amrex.v = 1 +amrex.verbose = 1 +amrex.regtest_reduction = 0 +amrex.signal_handling = 1 +amrex.throw_exception = 0 +amrex.call_addr2line = 1 +amrex.abort_on_unused_inputs = 0 +amrex.handle_sigsegv = 1 +amrex.handle_sigterm = 0 +amrex.handle_sigint = 1 +amrex.handle_sigabrt = 1 +amrex.handle_sigfpe = 1 +amrex.handle_sigill = 1 +amrex.fpe_trap_zero = 0 +amrex.fpe_trap_overflow = 0 +amrex.use_gpu_aware_mpi = 0 +amrex.the_arena_init_size = 0 +amrex.the_device_arena_init_size = 8388608 +amrex.the_managed_arena_init_size = 8388608 +amrex.the_pinned_arena_init_size = 8388608 +amrex.the_comms_arena_init_size = 8388608 +amrex.the_arena_release_threshold = 9223372036854775807 +amrex.the_device_arena_release_threshold = 9223372036854775807 +amrex.the_managed_arena_release_threshold = 9223372036854775807 +amrex.the_pinned_arena_release_threshold = 9223372036854775807 +amrex.the_comms_arena_release_threshold = 9223372036854775807 +amrex.the_async_arena_release_threshold = 9223372036854775807 +amrex.the_arena_is_managed = 0 +amrex.abort_on_out_of_gpu_memory = 0 +fab.init_snan = 1 +DistributionMapping.v = 0 +DistributionMapping.verbose = 0 +DistributionMapping.efficiency = 0.90000000000000002 +DistributionMapping.sfc_threshold = 0 +DistributionMapping.node_size = 0 +DistributionMapping.verbose_mapper = 0 +fab.initval = nan +fab.do_initval = 1 +fabarray.maxcomp = 25 +vismf.v = 0 +vismf.headerversion = 1 +vismf.groupsets = 0 +vismf.setbuf = 1 +vismf.usesingleread = 0 +vismf.usesinglewrite = 0 +vismf.checkfilepositions = 0 +vismf.usepersistentifstreams = 0 +vismf.usesynchronousreads = 0 +vismf.usedynamicsetselection = 1 +vismf.iobuffersize = 2097152 +vismf.allowsparsewrites = 1 +amrex.async_out = 0 +amrex.async_out_nfiles = 64 +amrex.vector_growth_factor = 1.5 +machine.verbose = 0 +machine.very_verbose = 0 +tiny_profiler.device_synchronize_around_region = 0 +tiny_profiler.verbose = 0 +tiny_profiler.v = 0 +amrex.use_profiler_syncs = 0 +geometry.coord_sys = 0 +amr.grid_eff = 0.69999999999999996 +amr.refine_grid_layout = 1 +amr.check_input = 1 diff --git a/Tests/test_files/Seamount/Seamount.i b/Tests/test_files/Seamount/Seamount.i new file mode 100644 index 00000000..f6dde426 --- /dev/null +++ b/Tests/test_files/Seamount/Seamount.i @@ -0,0 +1,71 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 10 +stop_time = 30000.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +# dims come from ROMS ana_grid, and must match vals in prob.cpp right now +geometry.prob_lo = 0. 0. -5000. +geometry.prob_hi = 320000. 320000. 0. + +amr.n_cell = 49 48 13 + +# periodic in x to match WRF setup +geometry.is_periodic = 1 1 0 +#ylo.type = "SlipWall" +#yhi.type = "SlipWall" +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +romsx.fixed_dt = 60.0 # Timestep size (seconds) +# NDTFAST = 30.0 # Number of baratropic steps => 300.0/30.0 = 10.0 +#romsx.fixed_fast_dt = 10.0 # Baratropic timestep size (seconds) +# romsx.fixed_fast_dt = 300.0 # Baratropic timestep size (seconds) testing value +romsx.fixed_ndtfast_ratio = 20 # Baratropic timestep size (seconds) + +romsx.flat_bathymetry=0 + +# DIAGNOSTICS & VERBOSITY +romsx.sum_interval = 1 # timesteps between computing mass +romsx.v = 0 # verbosity in ROMSX.cpp (0: none, 1: print boxes, etc, 2: print values) +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +romsx.check_file = chk # root name of checkpoint file +romsx.check_int = -57600 # number of timesteps between checkpoints + +# PLOTFILES +romsx.plot_file_1 = plt # prefix of plotfile name +romsx.plot_int_1 = 1 # number of timesteps between plotfiles +romsx.plot_vars_1 = omega salt temp x_velocity y_velocity z_velocity +romsx.plotfile_type = amrex + +# SOLVER CHOICE +romsx.use_coriolis = true +romsx.horizontal_advection_scheme = "upstream3" # upstream3 or centered4 +romsx.spatial_order = 2 + +# PROBLEM PARAMETERS (optional) +prob.R0 = 1027.0 +prob.S0 = 32.0 +prob.T0 = 10.0 + +# PROBLEM PARAMETERS (shear) +prob.rho_0 = 1.0 +prob.T_0 = 1.0 +prob.A_0 = 1.0 +prob.u_0 = 0.0 +prob.v_0 = 0.0 +prob.rad_0 = 0.25 +prob.z0 = 0.1 +prob.zRef = 80.0e-3 +prob.uRef = 8.0e-3 + +prob.prob_type = 12