Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

some cleanup but no change to answers #64

Merged
merged 1 commit into from
Dec 6, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 0 additions & 2 deletions Source/ROMSX.H
Original file line number Diff line number Diff line change
Expand Up @@ -357,7 +357,6 @@ public:
/** Advance the 3D variables */
void advance_3d (int lev,
amrex::MultiFab& mf_u , amrex::MultiFab& mf_v ,
amrex::MultiFab& mf_tempold , amrex::MultiFab& mf_saltold ,
amrex::MultiFab& mf_temp , amrex::MultiFab& mf_salt ,
std::unique_ptr<amrex::MultiFab>& mf_tempstore,
std::unique_ptr<amrex::MultiFab>& mf_saltstore,
Expand Down Expand Up @@ -582,7 +581,6 @@ public:
amrex::Array4<amrex::Real const> Dphi2,
amrex::Array4<amrex::Real> DC,
amrex::Array4<amrex::Real> FC,
amrex::Array4<amrex::Real> CF,
const int nnew);

void vert_mean_3d (const amrex::Box& bx,
Expand Down
44 changes: 20 additions & 24 deletions Source/TimeIntegration/ROMSX_advance_3d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ using namespace amrex;
void
ROMSX::advance_3d (int lev,
MultiFab& mf_u , MultiFab& mf_v ,
MultiFab& mf_tempold, MultiFab& mf_saltold,
MultiFab& mf_temp , MultiFab& mf_salt ,
std::unique_ptr<MultiFab>& mf_tempstore,
std::unique_ptr<MultiFab>& mf_saltstore,
Expand Down Expand Up @@ -77,7 +76,6 @@ ROMSX::advance_3d (int lev,
Array4<Real> const& Hvom = mf_Hvom->array(mfi);

Box bx = mfi.tilebox();
Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
Box gbx11 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,NGROW-1));
Box gbx21 = mfi.growntilebox(IntVect(NGROW,NGROW,NGROW-1));
Expand Down Expand Up @@ -160,53 +158,51 @@ ROMSX::advance_3d (int lev,
v(i,j,k) *= 2.0 / (Hz(i,j-1,k) + Hz(i,j,k));
});

{
amrex::Gpu::synchronize();
amrex::Gpu::LaunchSafeGuard lsg(true);

// NOTE: DC is only used as scratch in vert_visc_3d -- no need to pass or return a value
// NOTE: may not actually need to set these to zero

// Reset to zero on the box on which they'll be used
mf_DC[mfi].template setVal<RunOn::Device>(0.,gbx21);
fab_CF.template setVal<RunOn::Device>(0.,gbx21);
fab_CF.template setVal<RunOn::Device>(0.,gbx21);

#ifdef AMREX_USE_GPU
Gpu::synchronize();
#endif
vert_visc_3d(xbx,1,0,u,Hz,Hzk,oHz,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);

// Reset to zero
// Reset to zero on the box on which they'll be used
mf_DC[mfi].template setVal<RunOn::Device>(0.,gbx21);
fab_CF.template setVal<RunOn::Device>(0.,gbx21);
fab_CF.template setVal<RunOn::Device>(0.,gbx21);

vert_visc_3d(ybx,0,1,v,Hz,Hzk,oHz,AK,Akv,BC,DC,FC,CF,nnew,N,dt_lev);
}

// Reset to zero
mf_DC[mfi].template setVal<RunOn::Device>(0.,gbx21);
fab_CF.template setVal<RunOn::Device>(0.,gbx21);
// Reset to zero on the box on which they'll be used
mf_DC[mfi].template setVal<RunOn::Device>(0,xbx);
fab_CF.template setVal<RunOn::Device>(0.,xbx);

vert_mean_3d(xbx,1,0,u,Hz,DU_avg1,DC,CF,pm,nnew,N);

// Reset to zero
mf_DC[mfi].template setVal<RunOn::Device>(0.,gbx21);
fab_CF.template setVal<RunOn::Device>(0.,gbx21);
// Reset to zero on the box on which they'll be used
mf_DC[mfi].template setVal<RunOn::Device>(0.,ybx);
fab_CF.template setVal<RunOn::Device>(0.,ybx);

vert_mean_3d(ybx,0,1,v,Hz,DV_avg1,DC,CF,pn,nnew,N);

update_massflux_3d(gbx2,1,0,u,ubar,Huon,Hz,on_u,DU_avg1,DU_avg2,DC,FC,CF,nnew);
// Reset to zero on the box on which they'll be used
mf_DC[mfi].template setVal<RunOn::Device>(0.,grow(xbx,IntVect(0,0,1)));
fab_FC.template setVal<RunOn::Device>(0.,xbx);

update_massflux_3d(xbx,1,0,u,ubar,Huon,Hz,on_u,DU_avg1,DU_avg2,DC,FC,nnew);

update_massflux_3d(gbx2,0,1,v,vbar,Hvom,Hz,om_v,DV_avg1,DV_avg2,DC,FC,CF,nnew);
// Reset to zero on the box on which they'll be used
mf_DC[mfi].template setVal<RunOn::Device>(0.,grow(ybx,IntVect(0,0,1)));
fab_FC.template setVal<RunOn::Device>(0.,ybx);

update_massflux_3d(ybx,0,1,v,vbar,Hvom,Hz,om_v,DV_avg1,DV_avg2,DC,FC,nnew);
}

mf_Huon->FillBoundary(geom[lev].periodicity());
mf_Hvom->FillBoundary(geom[lev].periodicity());

for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
Array4<Real> const& tempold = (mf_tempold).array(mfi);
Array4<Real> const& saltold = (mf_saltold).array(mfi);

Array4<Real> const& temp = (mf_temp).array(mfi);
Array4<Real> const& salt = (mf_salt).array(mfi);

Expand Down
11 changes: 1 addition & 10 deletions Source/TimeIntegration/ROMSX_advance_3d_ml.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,6 @@ void ROMSX::advance_3d_ml (int lev, Real dt_lev)
#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

const BoxArray& ba = S_old.boxArray();
const DistributionMapping& dm = S_old.DistributionMap();
Expand All @@ -42,7 +36,7 @@ void ROMSX::advance_3d_ml (int lev, Real dt_lev)

const int ncomp = 1;

advance_3d(lev, mf_u, mf_v, mf_tempold, mf_saltold, mf_temp, mf_salt, vec_t3[lev], vec_s3[lev], vec_ru[lev], vec_rv[lev],
advance_3d(lev, mf_u, mf_v, mf_temp, mf_salt, vec_t3[lev], vec_s3[lev], vec_ru[lev], vec_rv[lev],
vec_DU_avg1[lev], vec_DU_avg2[lev],
vec_DV_avg1[lev], vec_DV_avg2[lev],
vec_ubar[lev], vec_vbar[lev],
Expand All @@ -59,9 +53,6 @@ void ROMSX::advance_3d_ml (int lev, 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());

Expand Down
83 changes: 31 additions & 52 deletions Source/TimeIntegration/ROMSX_update_massflux_3d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ using namespace amrex;
* @param[in ] Dphi_avg2
* @param[inout] DC
* @param[inout] FC
* @param[inout] CF
* @param[in ] nnew component of velocity
*/

Expand All @@ -26,83 +25,63 @@ ROMSX::update_massflux_3d (const Box& bx, const int ioff, const int joff,
Array4<Real const> Hz, Array4<Real> om_v_or_on_u,
Array4<Real const> Dphi_avg1,
Array4<Real const> Dphi_avg2, Array4<Real> DC,
Array4<Real> FC, Array4<Real> CF, const int nnew)
Array4<Real> FC, const int nnew)
{
const int Mn = Geom(0).Domain().size()[0];
const int Mm = Geom(0).Domain().size()[1];
// const int Mn = Geom(0).Domain().size()[0];
// const int Mm = Geom(0).Domain().size()[1];

auto N = Geom(0).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ

auto bxD=bx;
auto bx_g1z=bx;
bxD.makeSlab(2,0);
bx_g1z.grow(IntVect(0,0,1));

auto geomdata = Geom(0).data();
bool NSPeriodic = geomdata.isPeriodic(1);
bool EWPeriodic = geomdata.isPeriodic(0);
// auto geomdata = Geom(0).data();
// bool NSPeriodic = geomdata.isPeriodic(1);
// bool EWPeriodic = geomdata.isPeriodic(0);

//Copied depth of water column calculation from DepthStretchTransform
//Compute thicknesses of U-boxes DC(i,j,0:N-1), total depth of the water column DC(i,j,-1), and
// incorrect vertical mean CF(i,j,-1)

ParallelFor(bx_g1z, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
DC(i,j,k)=0.0;
CF(i,j,k)=0.0;
});

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
FC(i,j,k)=0.0;
});

Gpu::streamSynchronize();
//Compute thicknesses of U-boxes DC(i,j,0:N-1), total depth of the water column DC(i,j,-1)

//This takes advantage of Hz being an extra grow cell size
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
DC(i,j,k)=0.5*om_v_or_on_u(i,j,0)*(Hz(i,j,k)+Hz(i-ioff,j-joff,k));
});
{
DC(i,j,k)=0.5*om_v_or_on_u(i,j,0)*(Hz(i,j,k)+Hz(i-ioff,j-joff,k));
});

ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
for(int k=0; k<=N; k++) {
{
// Real CF = 0.;
for(int k=0; k<=N; k++) {
DC(i,j,-1)=DC(i,j,-1)+DC(i,j,k);
CF(i,j,-1)=CF(i,j,-1)+DC(i,j,k)*phi(i,j,k,nnew);
}
});
// CF += DC(i,j,k)*phi(i,j,k,nnew);
}

// Note this loop is in the opposite direction in k in ROMS but does not
// appear to affect results
DC(i,j,-1) = 1.0/DC(i,j,-1);
// CF = DC(i,j,-1) * (CF-Dphi_avg1(i,j,0));

ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
for (int k=0; k<=N; k++) {
if(k==0) {
DC(i,j,-1)=1.0/DC(i,j,-1);
CF(i,j,-1)=DC(i,j,-1)*(CF(i,j,-1)-Dphi_avg1(i,j,0));
}

if (!(NSPeriodic&&EWPeriodic))
{
if ( ( ((i<0)||(i>=Mn+1)) && !EWPeriodic ) || ( ((j<0)||(j>=Mm+1)) && !NSPeriodic ) )
{
phi(i,j,k) -= CF(i,j,-1);
}
}

// BOUNDARY CONDITIONS
// if (!(NSPeriodic&&EWPeriodic))
// {
// if ( ( ((i<0)||(i>=Mn+1)) && !EWPeriodic ) || ( ((j<0)||(j>=Mm+1)) && !NSPeriodic ) )
// {
// phi(i,j,k) -= CF;
// }
// }

Hphi(i,j,k) = 0.5 * (Hphi(i,j,k)+phi(i,j,k,nnew)*DC(i,j,k));
FC(i,j,0) = FC(i,j,0)+Hphi(i,j,k); //recursive
FC(i,j,0) += Hphi(i,j,k);
} // k
});

ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
{
FC(i,j,0) = DC(i,j,-1)*(FC(i,j,0)-Dphi_avg2(i,j,0)); //recursive
FC(i,j,0) = DC(i,j,-1)*(FC(i,j,0)-Dphi_avg2(i,j,0));
});

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Hphi(i,j,k) = Hphi(i,j,k)-DC(i,j,k)*FC(i,j,0);
Hphi(i,j,k) -= DC(i,j,k)*FC(i,j,0);
});

ParallelFor(bxD, [=] AMREX_GPU_DEVICE(int i, int j, int )
Expand Down