14 BL_PROFILE(
"REMORA::calc_stretch_coeffs()");
15 int nz = geom[0].Domain().length(2);
18 Real ds = 1.0_rt / Real(N);
20 Box bx(IntVect(0,0,0),IntVect(0,0,N));
28 auto s_w_dat =
s_w.data();
29 auto Cs_w_dat =
Cs_w.data();
30 auto s_r_dat =
s_r.data();
31 auto Cs_r_dat =
Cs_r.data();
33 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (
int ,
int ,
int k)
38 if (local_theta_s > 0.0_rt) {
39 Csur=(1.0_rt-std::cosh(local_theta_s*s_w_dat[k]))/
40 (std::cosh(local_theta_s)-1.0_rt);
42 Csur=-s_w_dat[k]*s_w_dat[k];
45 if (local_theta_b > 0.0_rt) {
46 Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
47 (1.0_rt-std::exp(-local_theta_b));
54 s_r_dat[k]=ds*(k-N+0.5_rt);
56 if (local_theta_s > 0.0_rt) {
57 Csur=(1.0_rt-std::cosh(local_theta_s*s_r_dat[k]))/
58 (std::cosh(local_theta_s)-1.0_rt);
60 Csur=-s_r_dat[k]*s_r_dat[k];
63 if (local_theta_b > 0.0_rt) {
64 Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
65 (1.0_rt-std::exp(-local_theta_b));
80 BL_PROFILE(
"REMORA::stretch_transform()");
81 std::unique_ptr<MultiFab>& mf_z_w =
vec_z_w[lev];
82 std::unique_ptr<MultiFab>& mf_z_r =
vec_z_r[lev];
83 std::unique_ptr<MultiFab>& mf_Hz =
vec_Hz[lev];
85 std::unique_ptr<MultiFab>& mf_Zt_avg1 =
vec_Zt_avg1[lev];
88 auto s_w_dat =
s_w.data();
89 auto Cs_w_dat =
Cs_w.data();
90 auto s_r_dat =
s_r.data();
91 auto Cs_r_dat =
Cs_r.data();
93 for ( MFIter mfi(*
cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
95 Array4<Real>
const& z_w = (mf_z_w)->array(mfi);
96 Array4<Real>
const& z_r = (mf_z_r)->array(mfi);
97 Array4<Real>
const& Hz = (mf_Hz)->array(mfi);
98 Array4<Real>
const& h = (mf_h)->array(mfi);
99 Array4<Real>
const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
100 Box bx = mfi.tilebox();
110 wgbx3.surroundingNodes(2);
112 const auto & geomdata = Geom(lev).data();
114 int nz = geom[lev].Domain().length(2);
119 Real ds = 1.0_rt / Real(N);
121 amrex::ParallelFor(wgbx3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
123 z_w(i,j,k) = h(i,j,0);
127 Gpu::streamSynchronize();
129 amrex::ParallelFor(wgbx3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
131 Real sc_r=ds*(k-N+0.5_rt);
135 Real cff1_r=Cs_r_dat[k];
138 Real cff1_w=Cs_w_dat[k];
140 Real hwater=h(i,j,0);
142 Real hinv=1.0_rt/(hc+hwater);
143 Real cff2_r=(cff_r+cff1_r*hwater)*hinv;
144 Real cff2_w=(cff_w+cff1_w*hwater)*hinv;
148 z_w(i,j,N)=Zt_avg1(i,j,0);
151 h(i,j,0,1) = Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
152 z_w(i,j,0) = h(i,j,0,1);
155 z_w(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
160 z_r(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_r;
164 Gpu::streamSynchronize();
166 amrex::ParallelFor(gbx3, [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
168 Hz(i,j,k)=z_w(i,j,k+1)-z_w(i,j,k);
172 vec_z_w[lev]->FillBoundary(geom[lev].periodicity());
173 vec_z_r[lev]->FillBoundary(geom[lev].periodicity());
174 vec_Hz[lev]->FillBoundary(geom[lev].periodicity());
177 for ( MFIter mfi(*
cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
179 Array4<Real>
const& z_w = (mf_z_w)->array(mfi);
180 Array4<Real>
const& z_phys_nd = (mf_z_phys_nd)->array(mfi);
182 Box z_w_box = Box(z_w);
183 auto const lo = amrex::lbound(z_w_box);
184 auto const hi = amrex::ubound(z_w_box);
194 Box bx = mfi.growntilebox(IntVect(
NGROW,
NGROW,0));
196 ParallelFor(convert(bx,IntVect(0,0,1)), [=] AMREX_GPU_DEVICE (
int i,
int j,
int k)
200 if ( i >= lo.x && i <= hi.x-1 && j >= lo.y && j <= hi.y-1 )
202 z_phys_nd(i,j,k)=0.25_rt*( z_w(i,j ,k) + z_w(i+1,j ,k) +
203 z_w(i,j+1,k) + z_w(i+1,j+1,k) );
205 int ii = std::min(std::max(i, lo.x), hi.x);
206 int jj = std::min(std::max(j, lo.y), hi.y);
207 z_phys_nd(i,j,k) = z_w(ii,jj,k);
213 ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (
int i,
int j,
int)
215 z_phys_nd(i,j,klo) = 2.0_rt * z_phys_nd(i,j,klo+1) - z_phys_nd(i,j,klo+2);