14 std::unique_ptr<MultiFab>& mf_z_w = vec_z_w[lev];
15 std::unique_ptr<MultiFab>& mf_z_r = vec_z_r[lev];
16 std::unique_ptr<MultiFab>& mf_s_r = vec_s_r[lev];
17 std::unique_ptr<MultiFab>& mf_Hz = vec_Hz[lev];
18 std::unique_ptr<MultiFab>& mf_h = vec_hOfTheConfusingName[lev];
19 std::unique_ptr<MultiFab>& mf_Zt_avg1 = vec_Zt_avg1[lev];
20 std::unique_ptr<MultiFab>& mf_z_phys_nd = vec_z_phys_nd[lev];
21 auto N_loop = Geom(lev).Domain().size()[2]-1;
23 for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
25 Array4<Real>
const& z_w = (mf_z_w)->array(mfi);
26 Array4<Real>
const& z_r = (mf_z_r)->array(mfi);
27 Array4<Real>
const& s_r = (mf_s_r)->array(mfi);
28 Array4<Real>
const& Hz = (mf_Hz)->array(mfi);
29 Array4<Real>
const& h = (mf_h)->array(mfi);
30 Array4<Real>
const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
31 Box bx = mfi.tilebox();
40 Box wgbx3 = gbx3.surroundingNodes(2);
42 const auto & geomdata = Geom(lev).data();
44 int nz = geom[lev].Domain().length(2);
48 Real hc=-min(geomdata.ProbHi(2),-solverChoice.tcline);
49 const auto local_theta_s = solverChoice.theta_s;
50 const auto local_theta_b = solverChoice.theta_b;
52 amrex::ParallelFor(wgbx3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
54 z_w(i,j,k) = h(i,j,0);
58 Gpu::streamSynchronize();
60 amrex::ParallelFor(gbx3D, [=] AMREX_GPU_DEVICE (
int i,
int j,
int )
62 for (
int k=-1; k<=N_loop; k++) {
80 Real cff_r, cff_w, cff1_r, cff1_w, cff2_r, cff2_w, Csur, Cbot;
81 Real sc_r,sc_w,Cs_r,Cs_w;
97 if (local_theta_s > 0.0_rt) {
98 Csur=(1.0_rt-std::cosh(local_theta_s*sc_w))/
99 (std::cosh(local_theta_s)-1.0_rt);
104 if (local_theta_b > 0.0_rt) {
105 Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
106 (1.0_rt-std::exp(-local_theta_b));
120 sc_r=ds*(k-N+0.5_rt);
122 if (local_theta_s > 0.0_rt) {
123 Csur=(1.0_rt-std::cosh(local_theta_s*sc_r))/
124 (std::cosh(local_theta_s)-1.0_rt);
129 if (local_theta_b > 0.0_rt) {
130 Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
131 (1.0_rt-std::exp(-local_theta_b));
138 if (i==0&&j==0&&k<N&&k>=0) {
146 Real hwater=h(i,j,0);
161 Real hinv=1.0_rt/(hc+hwater);
162 cff2_r=(cff_r+cff1_r*hwater)*hinv;
163 cff2_w=(cff_w+cff1_w*hwater)*hinv;
167 z_w(i,j,N-1)=Zt_avg1(i,j,0);
170 h(i,j,0,1) = Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
173 z_w(i,j,k-1)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
178 z_r(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_r;
183 Gpu::streamSynchronize();
185 amrex::ParallelFor(gbx3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
188 Hz(i,j,k)=z_w(i,j,k)+h(i,j,0);
190 Hz(i,j,k)=z_w(i,j,k)-z_w(i,j,k-1);
195 Real time = t_new[lev];
196 FillPatch(lev, time, *vec_z_w[lev], GetVecOfPtrs(vec_z_w));
197 FillPatch(lev, time, *vec_z_r[lev], GetVecOfPtrs(vec_z_r));
198 FillPatch(lev, time, *vec_s_r[lev], GetVecOfPtrs(vec_s_r));
199 FillPatch(lev, time, *vec_Hz[lev] , GetVecOfPtrs(vec_Hz));
202 for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
204 Array4<Real>
const& z_w = (mf_z_w)->array(mfi);
205 Array4<Real>
const& z_phys_nd = (mf_z_phys_nd)->array(mfi);
207 Box z_w_box = Box(z_w);
208 auto const lo = amrex::lbound(z_w_box);
209 auto const hi = amrex::ubound(z_w_box);
215 ParallelFor(Box(z_phys_nd), [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
219 int kk = (k == 0) ? hi.z : k-1;
220 if ( i >= lo.x && i <= hi.x-1 && j >= lo.y && j <= hi.y-1 )
222 z_phys_nd(i,j,k)=0.25*( z_w(i,j ,kk) + z_w(i+1,j ,kk) +
223 z_w(i,j+1,kk) + z_w(i+1,j+1,kk) );
225 int ii = std::min(std::max(i, lo.x), hi.x);
226 int jj = std::min(std::max(j, lo.y), hi.y);
227 z_phys_nd(i,j,k) = z_w(ii,jj,kk);
229 if (k == 0) z_phys_nd(i,j,k) *= -1.;
235 vec_z_phys_nd[lev]->FillBoundary(geom[lev].periodicity());
#define NGROW
Definition: IndexDefines.H:13
void stretch_transform(int lev)
Definition: DepthStretchTransform.H:12