Skip to content

Commit

Permalink
Merge pull request seahorce-scidac#72 from seahorce-scidac/more_to_ml
Browse files Browse the repository at this point in the history 
more cleanup towards multilevel
  • Loading branch information
@asalmgren
asalmgren authored Dec 10, 2023
2 parents c281b3e + 7fbe373 commit 5cefe81
Show file tree
Hide file tree
Showing 11 changed files with 153 additions and 193 deletions.
30 changes: 13 additions & 17 deletions Source/ROMSX.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::MultiFab>& mf_tempstore,
std::unique_ptr<amrex::MultiFab>& mf_saltstore,
std::unique_ptr<amrex::MultiFab>& mf_ru,
std::unique_ptr<amrex::MultiFab>& 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<amrex::MultiFab>& mf_DU_avg1,
std::unique_ptr<amrex::MultiFab>& mf_DU_avg2,
std::unique_ptr<amrex::MultiFab>& mf_DV_avg1,
std::unique_ptr<amrex::MultiFab>& mf_DV_avg2,
std::unique_ptr<amrex::MultiFab>& mf_ubar,
std::unique_ptr<amrex::MultiFab>& mf_vbar,
amrex::MultiFab& mf_AK , amrex::MultiFab& mf_DC,
amrex::MultiFab& mf_Hzk,
std::unique_ptr<amrex::MultiFab>& mf_Akv,
std::unique_ptr<amrex::MultiFab>& mf_Akt,
std::unique_ptr<amrex::MultiFab>& mf_Hz,
Expand Down Expand Up @@ -501,14 +497,14 @@ public:
int nrhs);

void rho_eos (const amrex::Box& bx,
amrex::Array4<amrex::Real> temp,
amrex::Array4<amrex::Real> salt,
amrex::Array4<amrex::Real> rho,
amrex::Array4<amrex::Real> rhoA,
amrex::Array4<amrex::Real> rhoS,
amrex::Array4<amrex::Real> Hz,
amrex::Array4<amrex::Real> z_w,
amrex::Array4<amrex::Real> h,
amrex::Array4<amrex::Real const> temp,
amrex::Array4<amrex::Real const> salt,
amrex::Array4<amrex::Real > rho,
amrex::Array4<amrex::Real > rhoA,
amrex::Array4<amrex::Real > rhoS,
amrex::Array4<amrex::Real const> Hz,
amrex::Array4<amrex::Real const> z_w,
amrex::Array4<amrex::Real const> h,
const int nrhs, const int N);

void prsgrd (const amrex::Box& bx,
Expand Down Expand Up @@ -566,7 +562,7 @@ public:
amrex::Array4<amrex::Real> v,
amrex::Array4<amrex::Real> Hv,
amrex::Array4<amrex::Real> om_v,
amrex::Array4<amrex::Real> Hz,
amrex::Array4<amrex::Real const> Hz,
const int nnew);

void update_massflux_3d (const amrex::Box& bx,
Expand Down
1 change: 0 additions & 1 deletion Source/ROMSX_SumIQ.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<const Real> vel_arr = cc_vel_mf.const_array(mfi);
const Array4<const Real> 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) +
Expand Down
14 changes: 11 additions & 3 deletions Source/TimeIntegration/ROMSX_Advance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
1 change: 0 additions & 1 deletion Source/TimeIntegration/ROMSX_TimeStepML.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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++)
Expand Down
92 changes: 47 additions & 45 deletions Source/TimeIntegration/ROMSX_advance_3d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<MultiFab>& mf_tempstore,
std::unique_ptr<MultiFab>& mf_saltstore,
std::unique_ptr<MultiFab>& mf_ru,
std::unique_ptr<MultiFab>& 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<MultiFab>& mf_DU_avg1,
std::unique_ptr<MultiFab>& mf_DU_avg2,
std::unique_ptr<MultiFab>& mf_DV_avg1,
std::unique_ptr<MultiFab>& mf_DV_avg2,
std::unique_ptr<MultiFab>& mf_ubar,
std::unique_ptr<MultiFab>& mf_vbar,
MultiFab& mf_AK, MultiFab& mf_DC,
MultiFab& mf_Hzk,
std::unique_ptr<MultiFab>& mf_Akv,
std::unique_ptr<MultiFab>& mf_Akt,
std::unique_ptr<MultiFab>& mf_Hz,
Expand All @@ -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<Real> const& u = mf_u.array(mfi);
Expand Down Expand Up @@ -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<Real> const& temp = (mf_temp).array(mfi);
Expand Down Expand Up @@ -283,8 +294,6 @@ ROMSX::advance_3d (int lev,
fab_on_u.template setVal<RunOn::Device>(dx[1],makeSlab(tbxp2,2,0));
fab_om_v.template setVal<RunOn::Device>(dx[0],makeSlab(tbxp2,2,0));

bool test_functionality=true;
if (test_functionality) {
//
//------------------------------------------------------------------------
// Vertically integrate horizontal mass flux divergence.
Expand All @@ -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());

Expand Down
55 changes: 19 additions & 36 deletions Source/TimeIntegration/ROMSX_advance_3d_ml.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

}
Loading

0 comments on commit 5cefe81

Please sign in to comment.