Skip to content

Commit

Permalink
Merge pull request #30 from jmsexton03/bathymetry_2d_boxes
Browse files Browse the repository at this point in the history 
Bathymetry 2d boxes
  • Loading branch information
@asalmgren
asalmgren authored Sep 20, 2023
2 parents 6880f89 + 571ee6c commit 9f4180e
Show file tree
Hide file tree
Showing 5 changed files with 219 additions and 24 deletions.
6 changes: 0 additions & 6 deletions Source/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,8 +66,6 @@ struct SolverChoice {

pp.query("flat_bathymetry", flat_bathymetry);

pp.query("test_vertical", test_vertical);

if(!flat_bathymetry) amrex::Warning("Flat bathymetry needs further testing");

// These default to true but are used for unit testing
Expand Down Expand Up @@ -198,9 +196,6 @@ struct SolverChoice {
//This accounts for the main 2d loop but may not include coupling and copying properly
pp.query("use_barotropic", use_barotropic);

//This accounts for the main 3d loop but may not include coupling and copying properly
pp.query("use_baroclinic", use_baroclinic);

static std::string ic_bc_type_string = "Ideal";
pp.query("ic_bc_type", ic_bc_type_string);

Expand Down Expand Up @@ -268,7 +263,6 @@ struct SolverChoice {
std::string pp_prefix {"romsx"};

bool flat_bathymetry = true;
bool test_vertical = true;

// Specify what additional physics/forcing modules we use
bool use_gravity = false;
Expand Down
5 changes: 3 additions & 2 deletions Source/TimeIntegration/ROMSX_Advance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -335,6 +335,9 @@ ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle
nstp, nrhs, N, dt_lev);
}


mf_W.FillBoundary(geom[lev].periodicity());

for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
Array4<Real> const& DC = mf_DC.array(mfi);
Expand Down Expand Up @@ -619,14 +622,12 @@ ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle
// because zeta may have changed
stretch_transform(lev);

if(solverChoice.use_baroclinic) {
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],
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_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());
Expand Down
207 changes: 205 additions & 2 deletions Source/TimeIntegration/ROMSX_advance_2d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,11 @@ ROMSX::advance_2d (int lev,
kstp-=1;
indx1-=1;
ptsk-=1;
auto ba = mf_h->boxArray();
auto dm = mf_h->DistributionMap();

MultiFab mf_DUon(ba,dm,1,IntVect(NGROW+1,NGROW+1,0));
MultiFab mf_DVom(ba,dm,1,IntVect(NGROW+1,NGROW+1,0));

for ( MFIter mfi(*mf_ru, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
Expand Down Expand Up @@ -173,8 +177,8 @@ ROMSX::advance_2d (int lev,
FArrayBox fab_gzeta2(tbxp2,1,The_Async_Arena());
FArrayBox fab_gzetaSA(tbxp2,1,The_Async_Arena());
FArrayBox fab_Dstp(tbxp2,1,The_Async_Arena());
FArrayBox fab_DUon(tbxp2,1,The_Async_Arena());
FArrayBox fab_DVom(tbxp2,1,The_Async_Arena());
FArrayBox & fab_DUon=mf_DUon[mfi];
FArrayBox & fab_DVom=mf_DVom[mfi];
FArrayBox fab_rhs_ubar(tbxp2,1,The_Async_Arena());
FArrayBox fab_rhs_vbar(tbxp2,1,The_Async_Arena());
FArrayBox fab_rhs_zeta(tbxp2uneven,1,The_Async_Arena());
Expand Down Expand Up @@ -281,6 +285,205 @@ ROMSX::advance_2d (int lev,
Real cff1=gbx2D.contains(i,j-1,0) ? cff*(Drhs(i,j,0)+Drhs(i,j-1,0)) : om_v(i,j,0)*Drhs(i,j,0);
DVom(i,j,0)=vbar(i,j,0,krhs)*cff1;
});
}
mf_DUon.FillBoundary(geom[lev].periodicity());
mf_DVom.FillBoundary(geom[lev].periodicity());

for ( MFIter mfi(*mf_ru, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{

Array4<Real> const& u = (mf_u).array(mfi);
Array4<Real> const& v = (mf_v).array(mfi);
Array4<Real> const& rhoS = (mf_rhoS).array(mfi);
Array4<Real> const& rhoA = (mf_rhoA).array(mfi);
Array4<Real> const& ubar = (mf_ubar)->array(mfi);
Array4<Real> const& vbar = (mf_vbar)->array(mfi);
Array4<Real> const& zeta = (mf_zeta)->array(mfi);
Array4<Real> const& h = (mf_h)->array(mfi);
Array4<Real> const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
Array4<Real> const& DU_avg1 = (mf_DU_avg1)->array(mfi);
Array4<Real> const& DU_avg2 = (mf_DU_avg2)->array(mfi);
Array4<Real> const& DV_avg1 = (mf_DV_avg1)->array(mfi);
Array4<Real> const& DV_avg2 = (mf_DV_avg2)->array(mfi);
Array4<Real> const& ru = (mf_ru)->array(mfi);
Array4<Real> const& rv = (mf_rv)->array(mfi);
Array4<Real> const& rufrc = (mf_rufrc)->array(mfi);
Array4<Real> const& rvfrc = (mf_rvfrc)->array(mfi);
Array4<Real> const& rubar = (mf_rubar)->array(mfi);
Array4<Real> const& rvbar = (mf_rvbar)->array(mfi);
Array4<Real> const& rzeta = (mf_rzeta)->array(mfi);
Array4<Real> const& visc2_p = (mf_visc2_p)->array(mfi);
Array4<Real> const& visc2_r = (mf_visc2_r)->array(mfi);

Box bx = mfi.tilebox();
Box gbx = mfi.growntilebox();
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 tbxp1 = bx;
Box tbxp11 = bx;
Box tbxp2 = bx;
tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
tbxp2.grow(IntVect(NGROW,NGROW,0));
tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
//make only gbx be grown to match multifabs
//gbx2.grow(IntVect(NGROW,NGROW,0));
//gbx1.grow(IntVect(NGROW-1,NGROW-1,0));
//gbx11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
Box ubxD = surroundingNodes(bx,0);
Box vbxD = surroundingNodes(bx,1);

Box bxD = bx;
bxD.makeSlab(2,0);
Box gbxD = gbx;
gbxD.makeSlab(2,0);
Box gbx1D = gbx1;
gbx1D.makeSlab(2,0);
Box gbx2D = gbx2;
gbx2D.makeSlab(2,0);

Box tbxp2D = tbxp2;
tbxp2D.makeSlab(2,0);


//bxD.makeSlab(2,0);
//ubxD.makeSlab(2,0);
//vbxD.makeSlab(2,0);
//Box gbx1D = bxD;
//Box gbx2D = bxD;
//gbx1D.grow(IntVect(NGROW-1,NGROW-1,0));
//gbx2D.grow(IntVect(NGROW,NGROW,0));

Box tbxp2uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-NGROW,bx.smallEnd(1)-NGROW,bx.smallEnd(2))),
IntVect(AMREX_D_DECL(bx.bigEnd(0)+NGROW-1,bx.bigEnd(1)+NGROW-1,bx.bigEnd(2))));

Box tbxp2unevenD = tbxp2uneven;
tbxp2unevenD.makeSlab(2,0);

BoxArray ba_gbx2uneven = intersect(BoxArray(tbxp2uneven), gbx);
AMREX_ASSERT((ba_gbx2uneven.size() == 1));
Box gbx2uneven = ba_gbx2uneven[0];
Box gbx2unevenD = gbx2uneven;
gbx2unevenD.makeSlab(2,0);
//AKA
//ubxD.setRange(2,0);
//vbxD.setRange(2,0);

FArrayBox fab_pn(tbxp2,1,The_Async_Arena());
FArrayBox fab_pm(tbxp2,1,The_Async_Arena());
FArrayBox fab_on_u(tbxp2,1,The_Async_Arena());
FArrayBox fab_om_v(tbxp2,1,The_Async_Arena());
FArrayBox fab_fomn(tbxp2,1,The_Async_Arena());
FArrayBox fab_Huon(tbxp2,1,The_Async_Arena()); //fab_Huon.setVal(0.0);
FArrayBox fab_Hvom(tbxp2,1,The_Async_Arena()); //fab_Hvom.setVal(0.0);
FArrayBox fab_oHz(tbxp11,1,The_Async_Arena()); //fab_oHz.setVal(0.0);
FArrayBox fab_om_u(tbxp2,1,amrex::The_Async_Arena());
FArrayBox fab_on_v(tbxp2,1,amrex::The_Async_Arena());
FArrayBox fab_om_r(tbxp2,1,amrex::The_Async_Arena());
FArrayBox fab_on_r(tbxp2,1,amrex::The_Async_Arena());
FArrayBox fab_om_p(tbxp2,1,amrex::The_Async_Arena());
FArrayBox fab_on_p(tbxp2,1,amrex::The_Async_Arena());

//step2d work arrays
FArrayBox fab_Drhs(tbxp2,1,The_Async_Arena());
FArrayBox fab_Drhs_p(tbxp2,1,The_Async_Arena());
FArrayBox fab_Dnew(tbxp2,1,The_Async_Arena());
FArrayBox fab_zwrk(tbxp2,1,The_Async_Arena());
FArrayBox fab_gzeta(tbxp2,1,The_Async_Arena());
FArrayBox fab_gzeta2(tbxp2,1,The_Async_Arena());
FArrayBox fab_gzetaSA(tbxp2,1,The_Async_Arena());
FArrayBox fab_Dstp(tbxp2,1,The_Async_Arena());
FArrayBox & fab_DUon=mf_DUon[mfi];
FArrayBox & fab_DVom=mf_DVom[mfi];
FArrayBox fab_rhs_ubar(tbxp2,1,The_Async_Arena());
FArrayBox fab_rhs_vbar(tbxp2,1,The_Async_Arena());
FArrayBox fab_rhs_zeta(tbxp2uneven,1,The_Async_Arena());
FArrayBox fab_zeta_new(tbxp2uneven,1,The_Async_Arena());

auto on_u=fab_on_u.array();
auto om_v=fab_om_v.array();
auto pn=fab_pn.array();
auto pm=fab_pm.array();
auto fomn=fab_fomn.array();
auto om_u=fab_om_u.array();
auto on_v=fab_on_v.array();
auto om_r=fab_om_r.array();
auto on_r=fab_on_r.array();
auto om_p=fab_om_p.array();
auto on_p=fab_on_p.array();

auto Drhs=fab_Drhs.array();
auto Drhs_p=fab_Drhs_p.array();
auto Dnew=fab_Dnew.array();
auto zwrk=fab_zwrk.array();
auto gzeta=fab_gzeta.array();
auto gzeta2=fab_gzeta2.array();
auto gzetaSA=fab_gzetaSA.array();
auto Dstp=fab_Dstp.array();
auto DUon=fab_DUon.array();
auto DVom=fab_DVom.array();
auto rhs_ubar=fab_rhs_ubar.array();
auto rhs_vbar=fab_rhs_vbar.array();
auto rhs_zeta=fab_rhs_zeta.array();
auto zeta_new=fab_zeta_new.array();

auto weight1 = vec_weight1.dataPtr();
auto weight2 = vec_weight2.dataPtr();
if ((verbose > 2) && predictor_2d_step && my_iif == 0) {
amrex::PrintToFile("ru_startadvance").SetPrecision(18)<<FArrayBox(ru)<<std::endl;
amrex::PrintToFile("rv_startadvance").SetPrecision(18)<<FArrayBox(rv)<<std::endl;
amrex::PrintToFile("u_startadvance2").SetPrecision(18)<<FArrayBox(u)<<std::endl;
amrex::PrintToFile("v_startadvance2").SetPrecision(18)<<FArrayBox(v)<<std::endl;
}

//From ana_grid.h and metrics.F
amrex::ParallelFor(tbxp2D,
[=] AMREX_GPU_DEVICE (int i, int j, int)
{
pm(i,j,0)=dxi[0];
pn(i,j,0)=dxi[1];
rhs_ubar(i,j,0)=0.0;
rhs_vbar(i,j,0)=0.0;
});

amrex::ParallelFor(tbxp2D,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
//Note: are the comment definitons right? Don't seem to match metrics.f90
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))
});
amrex::ParallelFor(tbxp2D,
[=] AMREX_GPU_DEVICE (int i, int j, int )
{

const auto prob_lo = geomdata.ProbLo();
const auto dx = geomdata.CellSize();

pm(i,j,0)=dxi[0];
pn(i,j,0)=dxi[1];
//defined UPWELLING
Real f0=-8.26e-5;
Real beta=0.0;
Real Esize=1000*(Mm);
Real y = prob_lo[1] + (j + 0.5) * dx[1];
Real f=f0+beta*(y-.5*Esize);
fomn(i,j,0)=f*(1.0/(pm(i,j,0)*pn(i,j,0)));
});

amrex::ParallelFor(tbxp2D,
[=] AMREX_GPU_DEVICE (int i, int j, int)
{
Drhs(i,j,0)=zeta(i,j,0,krhs)+h(i,j,0);
});
if(predictor_2d_step)
{
if(first_2d_step) {
Expand Down
16 changes: 5 additions & 11 deletions Source/TimeIntegration/ROMSX_advance_3d.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,6 @@ ROMSX::advance_3d (int lev,
} else {
cff=0.25*dt_lev*23.0/12.0;
}
if(solverChoice.test_vertical) {
amrex::ParallelFor(gbx2,
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Expand Down Expand Up @@ -248,12 +247,6 @@ ROMSX::advance_3d (int lev,
vbar(i,j,k,0) = DC(i,j,-1)*DV_avg1(i,j,0);
vbar(i,j,k,1) = vbar(i,j,k,0);
});
//amrex::ParallelFor(ubx,
//[=] AMREX_GPU_DEVICE (int i, int j, int k)
//{
// printf("%d %d %d %25.25g Huon after massflux\n", i, j, k, Huon(i,j,k));
// });
}
}
for ( MFIter mfi(mf_temp, TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
Expand Down Expand Up @@ -582,6 +575,7 @@ ROMSX::advance_3d (int lev,
}
Akt(i,j,k)=1e-6;
});

amrex::ParallelFor(amrex::makeSlab(tbxp2,2,0),
[=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Expand All @@ -593,10 +587,10 @@ ROMSX::advance_3d (int lev,
PrintToFile("temp_beforevvisc").SetPrecision(18) << FArrayBox(temp) << std::endl;
PrintToFile("salt_beforevvisc").SetPrecision(18) << FArrayBox(salt) << std::endl;
}
if(solverChoice.test_vertical) {
vert_visc_3d(gbx1,bx,0,0,temp,Hz,Hzk,oHz,AK,Akt,BC,DC,FC,CF,nnew,N,dt_lev);
vert_visc_3d(gbx1,bx,0,0,salt,Hz,Hzk,oHz,AK,Akt,BC,DC,FC,CF,nnew,N,dt_lev);
}

vert_visc_3d(gbx1,bx,0,0,temp,Hz,Hzk,oHz,AK,Akt,BC,DC,FC,CF,nnew,N,dt_lev);
vert_visc_3d(gbx1,bx,0,0,salt,Hz,Hzk,oHz,AK,Akt,BC,DC,FC,CF,nnew,N,dt_lev);

if (verbose > 2) {
PrintToFile("temp_aftervvisc").SetPrecision(18) << FArrayBox(temp) << std::endl;
PrintToFile("salt_aftervvisc").SetPrecision(18) << FArrayBox(salt) << std::endl;
Expand Down
9 changes: 6 additions & 3 deletions Source/Utils/DepthStretchTransform.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,14 +370,17 @@ ROMSX::stretch_transform (int lev)
cff2_w=(cff_w+cff1_w*hwater)*hinv;

if(k==0) {
z_w(i,j,N-1)=0.0;
// HACK: should actually be the normal expression with coeffs evaluated at k=N-1
z_w(i,j,N-1)=Zt_avg1(i,j,0);
}
else if(k==-1)
h(i,j,0,1)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
else
else {
z_w(i,j,k-1)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
if(k!=-1)
}
if(k!=-1) {
z_r(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_r;
}
}
}
});
Expand Down

0 comments on commit 9f4180e

Please sign in to comment.