diff --git a/Source/ROMSX.H b/Source/ROMSX.H index bb934936..b1ed19a9 100644 --- a/Source/ROMSX.H +++ b/Source/ROMSX.H @@ -357,20 +357,16 @@ public: /** Advance the 3D variables */ void advance_3d (int lev, - amrex::MultiFab& mf_u , amrex::MultiFab& mf_v , - amrex::MultiFab& mf_temp , amrex::MultiFab& mf_salt , - std::unique_ptr& mf_tempstore, - std::unique_ptr& mf_saltstore, - std::unique_ptr& mf_ru, - std::unique_ptr& mf_rv, + amrex::MultiFab& mf_u , amrex::MultiFab& mf_v, + amrex::MultiFab& mf_temp , amrex::MultiFab& mf_salt, + amrex::MultiFab* mf_tempstore, amrex::MultiFab* mf_saltstore, + amrex::MultiFab* mf_ru , amrex::MultiFab* mf_rv, std::unique_ptr& mf_DU_avg1, std::unique_ptr& mf_DU_avg2, std::unique_ptr& mf_DV_avg1, std::unique_ptr& mf_DV_avg2, std::unique_ptr& mf_ubar, std::unique_ptr& mf_vbar, - amrex::MultiFab& mf_AK , amrex::MultiFab& mf_DC, - amrex::MultiFab& mf_Hzk, std::unique_ptr& mf_Akv, std::unique_ptr& mf_Akt, std::unique_ptr& mf_Hz, @@ -501,14 +497,14 @@ public: int nrhs); void rho_eos (const amrex::Box& bx, - amrex::Array4 temp, - amrex::Array4 salt, - amrex::Array4 rho, - amrex::Array4 rhoA, - amrex::Array4 rhoS, - amrex::Array4 Hz, - amrex::Array4 z_w, - amrex::Array4 h, + amrex::Array4 temp, + amrex::Array4 salt, + 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); void prsgrd (const amrex::Box& bx, @@ -566,7 +562,7 @@ public: amrex::Array4 v, amrex::Array4 Hv, amrex::Array4 om_v, - amrex::Array4 Hz, + amrex::Array4 Hz, const int nnew); void update_massflux_3d (const amrex::Box& bx, diff --git a/Source/ROMSX_SumIQ.cpp b/Source/ROMSX_SumIQ.cpp index 7d0b5d85..5aabbdc8 100644 --- a/Source/ROMSX_SumIQ.cpp +++ b/Source/ROMSX_SumIQ.cpp @@ -30,7 +30,6 @@ ROMSX::sum_integrated_quantities(Real time) const Box& bx = mfi.tilebox(); const Array4< Real> kineng_arr = kineng_mf.array(mfi); const Array4 vel_arr = cc_vel_mf.const_array(mfi); - const Array4 cons_arr = cons_new[lev]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { kineng_arr(i,j,k) = 0.5 * ( vel_arr(i,j,k,0)*vel_arr(i,j,k,0) + vel_arr(i,j,k,1)*vel_arr(i,j,k,1) + diff --git a/Source/TimeIntegration/ROMSX_Advance.cpp b/Source/TimeIntegration/ROMSX_Advance.cpp index 3f67ebe1..375b795f 100644 --- a/Source/TimeIntegration/ROMSX_Advance.cpp +++ b/Source/TimeIntegration/ROMSX_Advance.cpp @@ -6,21 +6,29 @@ using namespace amrex; void ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle*/) { - BL_PROFILE("ROMSX::Advance()"); setup_step(lev, time, dt_lev); if (solverChoice.use_barotropic) { - bool predictor_2d_step=true; int nfast_counter=nfast + 1; + + //*************************************************** //Compute fast timestep from dt_lev and ratio + //*************************************************** Real dtfast_lev=dt_lev/Real(fixed_ndtfast_ratio); - for(int my_iif = 0; my_iif < nfast_counter; my_iif++) { + + //*************************************************** + //Advance nfast_counter steps of the 2d integrator + //*************************************************** + for (int my_iif = 0; my_iif < nfast_counter; my_iif++) { advance_2d_onestep(lev, dt_lev, dtfast_lev, my_iif, nfast_counter); } } + //*************************************************** + //Advance one step of the 3d integrator + //*************************************************** advance_3d_ml(lev, dt_lev); } diff --git a/Source/TimeIntegration/ROMSX_TimeStepML.cpp b/Source/TimeIntegration/ROMSX_TimeStepML.cpp index 76544028..68bbbbe3 100644 --- a/Source/TimeIntegration/ROMSX_TimeStepML.cpp +++ b/Source/TimeIntegration/ROMSX_TimeStepML.cpp @@ -65,7 +65,6 @@ ROMSX::timeStepML (Real time, int iteration) if (solverChoice.use_barotropic) { - bool predictor_2d_step=true; int nfast_counter=nfast + 1; for (int lev=0; lev <= finest_level; lev++) diff --git a/Source/TimeIntegration/ROMSX_advance_3d.cpp b/Source/TimeIntegration/ROMSX_advance_3d.cpp index d2e27e8f..ac0d2c0c 100644 --- a/Source/TimeIntegration/ROMSX_advance_3d.cpp +++ b/Source/TimeIntegration/ROMSX_advance_3d.cpp @@ -8,20 +8,16 @@ using namespace amrex; void ROMSX::advance_3d (int lev, - MultiFab& mf_u , MultiFab& mf_v , - MultiFab& mf_temp , MultiFab& mf_salt , - std::unique_ptr& mf_tempstore, - std::unique_ptr& mf_saltstore, - std::unique_ptr& mf_ru, - std::unique_ptr& mf_rv, + MultiFab& mf_u , MultiFab& mf_v , + MultiFab& mf_temp , MultiFab& mf_salt , + MultiFab* mf_tempstore, MultiFab* mf_saltstore, + MultiFab* mf_ru , MultiFab* mf_rv, std::unique_ptr& mf_DU_avg1, std::unique_ptr& mf_DU_avg2, std::unique_ptr& mf_DV_avg1, std::unique_ptr& mf_DV_avg2, std::unique_ptr& mf_ubar, std::unique_ptr& mf_vbar, - MultiFab& mf_AK, MultiFab& mf_DC, - MultiFab& mf_Hzk, std::unique_ptr& mf_Akv, std::unique_ptr& mf_Akt, std::unique_ptr& mf_Hz, @@ -47,6 +43,17 @@ ROMSX::advance_3d (int lev, // Because zeta may have changed stretch_transform(lev); + // These temporaries used to be made in advance_3d_ml and passed in; + // now we make them here + + const BoxArray& ba = mf_temp.boxArray(); + const DistributionMapping& dm = mf_temp.DistributionMap(); + + //Only used locally, probably should be rearranged into FArrayBox declaration + MultiFab mf_AK (ba,dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate + MultiFab mf_DC (ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate + MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate + for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& u = mf_u.array(mfi); @@ -214,6 +221,10 @@ ROMSX::advance_3d (int lev, mf_Huon->FillBoundary(geom[lev].periodicity()); mf_Hvom->FillBoundary(geom[lev].periodicity()); + // ************************************************************************ + // This should fill both temp and salt with temp/salt currently in cons_old + // ************************************************************************ + for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { Array4 const& temp = (mf_temp).array(mfi); @@ -283,8 +294,6 @@ ROMSX::advance_3d (int lev, fab_on_u.template setVal(dx[1],makeSlab(tbxp2,2,0)); fab_om_v.template setVal(dx[0],makeSlab(tbxp2,2,0)); - bool test_functionality=true; - if (test_functionality) { // //------------------------------------------------------------------------ // Vertically integrate horizontal mass flux divergence. @@ -295,56 +304,49 @@ ROMSX::advance_3d (int lev, Box gbx1D = gbx1; gbx1D.makeSlab(2,0); - ParallelFor(gbx1D, - [=] AMREX_GPU_DEVICE (int i, int j, int ) - { // Starting with zero vertical velocity at the bottom, integrate // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng)) // contains the vertical velocity at the free-surface, d(zeta)/d(t). // Notice that barotropic mass flux divergence is not used directly. - // - //W(i,j,-1)=0.0; - int k=0; - W(i,j,k) = - (Huon(i+1,j,k)-Huon(i,j,k)) - (Hvom(i,j+1,k)-Hvom(i,j,k)); - for(k=1;k<=N;k++) { - W(i,j,k) = W(i,j,k-1) - (Huon(i+1,j,k)-Huon(i,j,k)) - (Hvom(i,j+1,k)-Hvom(i,j,k)); - } - }); + ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int ) + { + W(i,j,0) = - (Huon(i+1,j,0)-Huon(i,j,0)) - (Hvom(i,j+1,0)-Hvom(i,j,0)); + + for (int k=1; k<=N; k++) { + W(i,j,k) = W(i,j,k-1) - (Huon(i+1,j,k)-Huon(i,j,k)) - (Hvom(i,j+1,k)-Hvom(i,j,k)); + } + }); - ParallelFor(gbx1, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { // Starting with zero vertical velocity at the bottom, integrate // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng)) // contains the vertical velocity at the free-surface, d(zeta)/d(t). // Notice that barotropic mass flux divergence is not used directly. // - Real wrk_i=W(i,j,N)/(z_w(i,j,N)+h(i,j,0,0)); + ParallelFor(gbx1, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real wrk_ij = W(i,j,N) / (z_w(i,j,N)+h(i,j,0,0)); - if(k!=N) { - W(i,j,k) = W(i,j,k)- wrk_i*(z_w(i,j,k)+h(i,j,0,0)); - } - }); + if(k!=N) { + W(i,j,k) -= wrk_ij * (z_w(i,j,k)+h(i,j,0,0)); + } + }); - // probably not the most efficient way - ParallelFor(gbx1, - [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (k == N) { + ParallelFor(makeSlab(gbx1,2,N), + [=] AMREX_GPU_DEVICE (int i, int j, int) + { W(i,j,N) = 0.0; - } - }); - - // - //----------------------------------------------------------------------- - // rhs_t_3d - //----------------------------------------------------------------------- - // - rhs_t_3d(bx, gbx, temp, tempstore, Huon, Hvom, Hz, oHz, pn, pm, W, FC, nrhs, nnew, N,dt_lev); - rhs_t_3d(bx, gbx, salt, saltstore, Huon, Hvom, Hz, oHz, pn, pm, W, FC, nrhs, nnew, N,dt_lev); - } + }); + + // + //----------------------------------------------------------------------- + // rhs_t_3d + //----------------------------------------------------------------------- + // + rhs_t_3d(bx, gbx, temp, tempstore, Huon, Hvom, Hz, oHz, pn, pm, W, FC, nrhs, nnew, N,dt_lev); + rhs_t_3d(bx, gbx, salt, saltstore, Huon, Hvom, Hz, oHz, pn, pm, W, FC, nrhs, nnew, N,dt_lev); } // mfi + mf_temp.FillBoundary(geom[lev].periodicity()); mf_salt.FillBoundary(geom[lev].periodicity()); diff --git a/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp b/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp index 2ebb3f1e..f3edf691 100644 --- a/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp +++ b/Source/TimeIntegration/ROMSX_advance_3d_ml.cpp @@ -2,72 +2,55 @@ using namespace amrex; -// advance a single 3D level for a single time step +// Advance a single 3D level for a single time step void ROMSX::advance_3d_ml (int lev, Real dt_lev) { - MultiFab& S_old = *cons_old[lev]; - MultiFab& S_new = *cons_new[lev]; - - MultiFab& U_old = *xvel_old[lev]; - MultiFab& V_old = *yvel_old[lev]; - - MultiFab& U_new = *xvel_new[lev]; - MultiFab& V_new = *yvel_new[lev]; - - MultiFab mf_u(U_new, amrex::make_alias, 0, 1); - MultiFab mf_v(V_new, amrex::make_alias, 0, 1); + // Fill in three ways: 1) interpolate from coarse grid if lev > 0; 2) fill from physical boundaries; + // 3) fine-fine fill of ghost cells with FillBoundary call + FillPatch(lev, t_old[lev], cons_old[lev], cons_old); + FillPatch(lev, t_old[lev], xvel_old[lev], xvel_old); + FillPatch(lev, t_old[lev], yvel_old[lev], yvel_old); + FillPatch(lev, t_old[lev], zvel_old[lev], zvel_old); - MultiFab mf_temp(S_new, amrex::make_alias, Temp_comp, 1); + MultiFab mf_temp(*cons_new[lev], amrex::make_alias, Temp_comp, 1); #ifdef ROMSX_USE_SALINITY - MultiFab mf_salt(S_new, amrex::make_alias, Salt_comp, 1); + MultiFab mf_salt(*cons_new[lev], amrex::make_alias, Salt_comp, 1); #else - MultiFab mf_salt(S_new, amrex::make_alias, Temp_comp, 1); + MultiFab mf_salt(*cons_new[lev], amrex::make_alias, Temp_comp, 1); #endif - const BoxArray& ba = S_old.boxArray(); - const DistributionMapping& dm = S_old.DistributionMap(); - - //Only used locally, probably should be rearranged into FArrayBox declaration - MultiFab mf_AK(ba,dm,1,IntVect(NGROW,NGROW,0)); //2d missing j coordinate - MultiFab mf_DC(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate - MultiFab mf_Hzk(ba,dm,1,IntVect(NGROW,NGROW,NGROW-1)); //2d missing j coordinate - auto N = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ const int ncomp = 1; - advance_3d(lev, mf_u, mf_v, mf_temp, mf_salt, vec_t3[lev], vec_s3[lev], vec_ru[lev], vec_rv[lev], + advance_3d(lev, *xvel_new[lev], *yvel_new[lev], + mf_temp, mf_salt, + vec_t3[lev].get(), vec_s3[lev].get(), + vec_ru[lev].get(), vec_rv[lev].get(), vec_DU_avg1[lev], vec_DU_avg2[lev], vec_DV_avg1[lev], vec_DV_avg2[lev], vec_ubar[lev], vec_vbar[lev], - mf_AK, mf_DC, - mf_Hzk, vec_Akv[lev], vec_Akt[lev], vec_Hz[lev], vec_Huon[lev], vec_Hvom[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); - U_new.FillBoundary(geom[lev].periodicity()); - V_new.FillBoundary(geom[lev].periodicity()); - - U_old.FillBoundary(geom[lev].periodicity()); - V_old.FillBoundary(geom[lev].periodicity()); - vec_ubar[lev]->FillBoundary(geom[lev].periodicity()); vec_vbar[lev]->FillBoundary(geom[lev].periodicity()); - mf_temp.FillBoundary(geom[lev].periodicity()); - mf_salt.FillBoundary(geom[lev].periodicity()); - vec_t3[lev]->FillBoundary(geom[lev].periodicity()); vec_s3[lev]->FillBoundary(geom[lev].periodicity()); // Fill in three ways: 1) interpolate from coarse grid if lev > 0; 2) fill from physical boundaries; // 3) fine-fine fill of ghost cells with FillBoundary call + // Note that we need the fine-fine and physical bc's in order to correctly move the particles FillPatch(lev, t_new[lev], cons_new[lev], cons_new); FillPatch(lev, t_new[lev], xvel_new[lev], xvel_new); FillPatch(lev, t_new[lev], yvel_new[lev], yvel_new); FillPatch(lev, t_new[lev], zvel_new[lev], zvel_new); #ifdef ROMSX_USE_PARTICLES + //*************************************************** + //Advance particles + //*************************************************** particleData.advance_particles(lev, dt_lev, xvel_new[lev], yvel_new[lev], zvel_new[lev], vec_z_phys_nd); #endif - } diff --git a/Source/TimeIntegration/ROMSX_rho_eos.cpp b/Source/TimeIntegration/ROMSX_rho_eos.cpp index 67ffdc68..52eda244 100644 --- a/Source/TimeIntegration/ROMSX_rho_eos.cpp +++ b/Source/TimeIntegration/ROMSX_rho_eos.cpp @@ -19,19 +19,18 @@ using namespace amrex; void ROMSX::rho_eos (const Box& bx, - Array4 temp , Array4 salt, - Array4 rho, - Array4 rhoA, - Array4 rhoS, - Array4 Hz, - Array4 z_w, - Array4 h, + Array4 temp, + Array4 salt, + Array4 rho, + Array4 rhoA, + Array4 rhoS, + Array4 Hz, + Array4 z_w, + Array4 h, const int nrhs, const int N) { - auto bxD=bx; - bxD.makeSlab(2,0); - //hardcode these for now instead of reading them from inputs + // Hardcode these for now instead of reading them from inputs Real T0=14.0; Real S0=35.0; Real R0=1027; @@ -39,6 +38,7 @@ ROMSX::rho_eos (const Box& bx, Real Scoef=0.0; Real rho0=1025.0; + AMREX_ASSERT(bx.smallEnd(2) == 0 && bx.bigEnd(2) == N); // //======================================================================= // Linear equation of state. @@ -52,12 +52,8 @@ ROMSX::rho_eos (const Box& bx, amrex::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)=rho(i,j,k)+ - R0*Scoef*(salt(i,j,k,nrhs)-S0); - rho(i,j,k)=rho(i,j,k)-1000.0_rt; -// pden(i,j,k)=rho(i,j,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; }); // @@ -66,43 +62,24 @@ ROMSX::rho_eos (const Box& bx, // used (rhoS) in barotropic pressure gradient. //----------------------------------------------------------------------- // - amrex::ParallelFor(bxD, - [=] AMREX_GPU_DEVICE (int i, int j, int ) - { - //printf("%d %d %15.15g Hzstart\n", i,j, Hz(i,j,N)); - Real cff1=rho(i,j,N)*Hz(i,j,N); - rhoS(i,j,0)=0.5_rt*cff1*Hz(i,j,N); - rhoA(i,j,0)=cff1; - //printf("%d %d %d %15.15g cff1\n", i,j,N, cff1); - //printf("%d %d %d %15.15g rhoN\n", i,j,N, rho(i,j,N)); - //printf("%d %d %d %15.15g rhoA\n", i,j,N, rhoA(i,j,0)); - //printf("%d %d %d %15.15g rhoS\n", i,j,N, rhoS(i,j,0)); - //printf("%d %d %d %15.15g Hzend \n", i,j,N, Hz(i,j,N)); - //printf("%d %d %15.15g Hzend\n", i,j, Hz(i,j,N)); - //printf("%d %d %15.15g %15.15g %15.15g %15.15g %15.15g cff1 rhoN rhoA rhoS Hz rhoeos\n", - // i,j, cff1, rho(i,j,N), rhoA(i,j,0), rhoS(i,j,0), Hz(i,j,N)); - }); - AMREX_ASSERT(bx.smallEnd(2) == 0 && - bx.bigEnd(2) == N); - amrex::ParallelFor(bxD, - [=] AMREX_GPU_DEVICE (int i, int j, int) - { - for (int k = 1; k <= N; ++k) { - Real cff1=rho(i,j,N-k)*Hz(i,j,N-k); - rhoS(i,j,0)=rhoS(i,j,0)+Hz(i,j,N-k)*(rhoA(i,j,0)+0.5_rt*cff1); - rhoA(i,j,0)=rhoA(i,j,0)+cff1; - //printf("%d %d %d %15.15g %15.15g %15.15g %15.15g cff1 rhoA rhoS Hz rhoeos\n", i,j,k, cff1, rhoA(i,j,0), rhoS(i,j,0), Hz(i,j,N-k)); - } - }); - amrex::ParallelFor(bxD, - [=] AMREX_GPU_DEVICE (int i, int j, int ) + Real cff2 =1.0_rt/rho0; + + amrex::ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int ) { - Real cff2=1.0_rt/rho0; - Real cff1=1.0_rt/(z_w(i,j,N)+h(i,j,0,0)); - rhoA(i,j,0)=cff2*cff1*rhoA(i,j,0); - rhoS(i,j,0)=2.0_rt*cff1*cff1*cff2*rhoS(i,j,0); - //printf("%d %d %15.15g %15.15g %15.15g %15.15g z_wN z_w0\n", i,j, z_w(i,j,N), z_w(i,j,0), h(i,j,0,0),h(i,j,0,1)); - //printf("%d %d %15.15g %15.15g %15.15g %15.15g cff2 cff1 rhoA rhoS rhoeos\n", i,j,cff2, cff1, rhoA(i,j,0), rhoS(i,j,0)); - }); + Real cff0 = rho(i,j,N)*Hz(i,j,N); + rhoS(i,j,0) = 0.5_rt*cff0*Hz(i,j,N); + rhoA(i,j,0) = cff0; + + for (int k = 1; k <= N; ++k) { + Real cff1=rho(i,j,N-k)*Hz(i,j,N-k); + rhoS(i,j,0) += Hz(i,j,N-k)*(rhoA(i,j,0)+0.5_rt*cff1); + rhoA(i,j,0) += cff1; + } + + Real cff11 =1.0_rt/(z_w(i,j,N)+h(i,j,0,0)); + rhoA(i,j,0) *= cff2*cff11; + + rhoS(i,j,0) *= 2.0_rt*cff11*cff11*cff2; + }); } diff --git a/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp b/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp index 6418ea6d..954f20e7 100644 --- a/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp +++ b/Source/TimeIntegration/ROMSX_rhs_uv_2d.cpp @@ -139,7 +139,7 @@ ROMSX::rhs_uv_2d (const Box& xbx, const Box& ybx, }); amrex::ParallelFor(growLo(ybx,1,1), - [=] AMREX_GPU_DEVICE (int i, int j, int k) + [=] 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); Real vee_jp1 = vbar(i,j ,0,krhs)-2.0*vbar(i,j+1,0,krhs)+vbar(i,j+2,0,krhs); diff --git a/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp b/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp index ff8f73f4..a3ec8f55 100644 --- a/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp +++ b/Source/TimeIntegration/ROMSX_set_massflux_3d.cpp @@ -17,7 +17,7 @@ using namespace amrex; 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, const int nrhs) { // //----------------------------------------------------------------------- diff --git a/Source/TimeIntegration/ROMSX_set_weights.cpp b/Source/TimeIntegration/ROMSX_set_weights.cpp index 528dcb1f..8e3a22b1 100644 --- a/Source/TimeIntegration/ROMSX_set_weights.cpp +++ b/Source/TimeIntegration/ROMSX_set_weights.cpp @@ -10,7 +10,6 @@ using namespace amrex; void ROMSX::set_weights (int /*lev*/) { Real gamma, scale; - Real cff1, cff2; Real wsum, shift, cff; //HACK should possibly store fixed_ndtfast elsewhere @@ -165,6 +164,8 @@ void ROMSX::set_weights (int /*lev*/) { // // Report weights. // +#if 0 + Real cff1, cff2; if (ParallelDescriptor::IOProcessor()) { Print().SetPrecision(18)<nComp()); - FillPatchNoPhysBC(lev, time, U_old, {xvel_old[lev]}, 0, 1); - FillPatchNoPhysBC(lev, time, V_old, {yvel_old[lev]}, 0, 1); - FillPatchNoPhysBC(lev, time, W_old, {zvel_old[lev]}, 0, 1); + int nvars = S_old.nComp(); + + // Fill ghost cells/faces at old time + FillPatch(lev, time, cons_old[lev], cons_old); + FillPatch(lev, time, xvel_old[lev], xvel_old); + FillPatch(lev, time, yvel_old[lev], yvel_old); + FillPatch(lev, time, zvel_old[lev], zvel_old); ////////// //pre_step3d corrections to boundaries const BoxArray& ba = S_old.boxArray(); const DistributionMapping& dm = S_old.DistributionMap(); - int nvars = S_old.nComp(); - const int ncomp = 1; const int nrhs = ncomp-1; const int nnew = ncomp-1; @@ -54,11 +55,7 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) std::unique_ptr& mf_h = vec_hOfTheConfusingName[lev]; //Consider passing these into the advance function or renaming relevant things - /* - MultiFab mf_u(ba,dm,1,IntVect(NGROW,NGROW,0)); - MultiFab mf_v(ba,dm,1,IntVect(NGROW,NGROW,0)); - MultiFab mf_uold(ba,dm,1,IntVect(NGROW,NGROW,0)); - MultiFab mf_vold(ba,dm,1,IntVect(NGROW,NGROW,0));*/ + MultiFab mf_u(U_new, amrex::make_alias, 0, 1); MultiFab mf_v(V_new, amrex::make_alias, 0, 1); MultiFab mf_uold(U_old, amrex::make_alias, 0, 1); @@ -98,14 +95,10 @@ ROMSX::setup_step (int lev, Real time, Real dt_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 - /* - mf_u.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); - mf_v.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); - mf_uold.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); - mf_vold.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0)));*/ mf_rho.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); mf_rhoS->setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); mf_rhoA->setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); + mf_w.setVal(0); mf_DC.setVal(0); mf_w.setVal(0.e34,IntVect(AMREX_D_DECL(NGROW-1,NGROW-1,0))); @@ -155,21 +148,21 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) //MFIter::allowMultipleMFIters(true); for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - Array4 const& h = (vec_hOfTheConfusingName[lev])->array(mfi); - Array4 const& Hz = (vec_Hz[lev])->array(mfi); - Array4 const& Huon = (vec_Huon[lev])->array(mfi); - Array4 const& Hvom = (vec_Hvom[lev])->array(mfi); - Array4 const& z_w = (mf_z_w)->array(mfi); - Array4 const& uold = (mf_uold).array(mfi); - Array4 const& vold = (mf_vold).array(mfi); - 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).array(mfi); - Array4 const& saltold = (mf_saltold).array(mfi); - Array4 const& rdrag = (mf_rdrag)->array(mfi); - Array4 const& bustr = (mf_bustr)->array(mfi); - Array4 const& bvstr = (mf_bvstr)->array(mfi); + Array4 const& h = (vec_hOfTheConfusingName[lev])->const_array(mfi); + Array4 const& Hz = (vec_Hz[lev])->const_array(mfi); + Array4 const& Huon = (vec_Huon[lev])->array(mfi); + Array4 const& Hvom = (vec_Hvom[lev])->array(mfi); + Array4 const& z_w = (mf_z_w)->const_array(mfi); + Array4 const& uold = (mf_uold).array(mfi); + Array4 const& vold = (mf_vold).array(mfi); + 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); Box bx = mfi.tilebox(); Box ubx = Box(uold); @@ -256,8 +249,7 @@ ROMSX::setup_step (int lev, Real time, Real dt_lev) 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 ) + amrex::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) @@ -500,11 +492,13 @@ ROMSX::setup_step (int lev, Real time, Real dt_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()); + vec_t3[lev]->FillBoundary(geom[lev].periodicity()); vec_s3[lev]->FillBoundary(geom[lev].periodicity()); + vec_Huon[lev]->FillBoundary(geom[lev].periodicity()); vec_Hvom[lev]->FillBoundary(geom[lev].periodicity()); - }