REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_gls.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] lev level to operate on
7 * @param[inout] mf_gls turbulent generic length scale
8 * @param[inout] mf_tke turbulent kinetic energy
9 * @param[in ] mf_W vertical velocity
10 * @param[in ] mf_msku land-sea mask on u points
11 * @param[in ] mf_mskv land-sea mask on v points
12 * @param[in ] nstp index of last time step in gls and tke MultiFabs
13 * @param[in ] nnew index of time step to update in gls and tke MultiFabs
14 * @param[in ] iic which time step we're on
15 * @param[in ] ntfirst what is the first time step?
16 * @param[in ] N number of vertical levels
17 * @param[in ] dt_lev time step at this level
18 */
19void
20REMORA::gls_prestep (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
21 MultiFab& mf_W, MultiFab* mf_msku, MultiFab* mf_mskv,
22 const int nstp, const int nnew,
23 const int iic, const int ntfirst, const int N, const Real dt_lev)
24{
25 BL_PROFILE("REMORA::gls_prestep()");
26 // temps: grad, gradL, XF, FX, FXL, EF, FE, FEL
27 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
28 Array4<Real> const& gls = mf_gls->array(mfi);
29 Array4<Real> const& tke = mf_tke->array(mfi);
30 Array4<Real const> const& W = mf_W.const_array(mfi);
31
32 Array4<Real const> const& Huon = vec_Huon[lev]->const_array(mfi);
33 Array4<Real const> const& Hvom = vec_Hvom[lev]->const_array(mfi);
34 Array4<Real const> const& Hz = vec_Hz[lev]->const_array(mfi);
35 Array4<Real const> const& pm = vec_pm[lev]->const_array(mfi);
36 Array4<Real const> const& pn = vec_pn[lev]->const_array(mfi);
37 Array4<Real const> const& msku = mf_msku->const_array(mfi);
38 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
39
40 Box bx = mfi.tilebox();
41 Box xbx = surroundingNodes(bx,0);
42 Box ybx = surroundingNodes(bx,1);
43
44 Box xbx_hi = growHi(xbx,0,1);
45
46 Box ybx_hi = growHi(ybx,0,1);
47
48 const Box& domain = geom[lev].Domain();
49 const auto dlo = amrex::lbound(domain);
50 const auto dhi = amrex::ubound(domain);
51
52 GeometryData const& geomdata = geom[0].data();
53 bool is_periodic_in_x = geomdata.isPeriodic(0);
54 bool is_periodic_in_y = geomdata.isPeriodic(1);
55
56 int ncomp = 1;
57 Vector<BCRec> bcrs_x(ncomp);
58 Vector<BCRec> bcrs_y(ncomp);
59 amrex::setBC(xbx,domain,xvel_bc(),0,1,domain_bcs_type,bcrs_x);
60 amrex::setBC(ybx,domain,yvel_bc(),0,1,domain_bcs_type,bcrs_y);
61
62 FArrayBox fab_XF(xbx_hi, 1, amrex::The_Async_Arena()); fab_XF.template setVal<RunOn::Device>(zero);
63 FArrayBox fab_FX(xbx_hi, 1, amrex::The_Async_Arena()); fab_FX.template setVal<RunOn::Device>(zero);
64 FArrayBox fab_FXL(xbx_hi, 1, amrex::The_Async_Arena()); fab_FXL.template setVal<RunOn::Device>(zero);
65 FArrayBox fab_EF(ybx_hi, 1, amrex::The_Async_Arena()); fab_EF.template setVal<RunOn::Device>(zero);
66 FArrayBox fab_FE(ybx_hi, 1, amrex::The_Async_Arena()); fab_FE.template setVal<RunOn::Device>(zero);
67 FArrayBox fab_FEL(ybx_hi, 1, amrex::The_Async_Arena()); fab_FEL.template setVal<RunOn::Device>(zero);
68 FArrayBox fab_Hz_half(bx, 1, amrex::The_Async_Arena()); fab_Hz_half.template setVal<RunOn::Device>(zero);
69 FArrayBox fab_CF(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_CF.template setVal<RunOn::Device>(zero);
70 FArrayBox fab_FC(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_FC.template setVal<RunOn::Device>(zero);
71 FArrayBox fab_FCL(convert(bx,IntVect(0,0,0)), 1, amrex::The_Async_Arena()); fab_FCL.template setVal<RunOn::Device>(zero);
72
73 auto XF = fab_XF.array();
74 auto FX = fab_FX.array();
75 auto FXL = fab_FXL.array();
76 auto EF = fab_EF.array();
77 auto FE = fab_FE.array();
78 auto FEL = fab_FEL.array();
79 auto Hz_half = fab_Hz_half.array();
80 auto CF = fab_CF.array();
81 auto FC = fab_FC.array();
82 auto FCL = fab_FCL.array();
83
84 // need XF/FX/FXL from [xlo to xhi] by [ylo to yhi ] on u points
85 ParallelFor(grow(xbx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
86 {
87 Real grad_im1 = (tke(i-1,j,k,nstp) - tke(i-2,j,k,nstp)) * msku(i-1,j,0);
88 Real grad_ip1 = (tke(i+1,j,k,nstp) - tke(i ,j,k,nstp)) * msku(i+1,j,0);
89
90 Real gradL_im1 = (gls(i-1,j,k,nstp) - gls(i-2,j,k,nstp)) * msku(i-1,j,0);
91 Real gradL_ip1 = (gls(i+1,j,k,nstp) - gls(i ,j,k,nstp)) * msku(i+1,j,0);
92
93 // Adjust boundaries
94 // TODO: Make sure indices match with what ROMS does
95 if (i == dlo.x-1 && !is_periodic_in_x) {
96 grad_im1 = tke(i,j,k,nstp) - tke(i-1,j,k,nstp);
97 gradL_im1 = gls(i,j,k,nstp) - gls(i-1,j,k,nstp);
98 }
99 else if (i == dhi.x+1 && !is_periodic_in_x) {
100 grad_ip1 = tke(i,j,k,nstp) - tke(i-1,j,k,nstp);
101 gradL_ip1 = gls(i,j,k,nstp) - gls(i-1,j,k,nstp);
102 }
103 Real cff = one/Real(6.0);
104 XF(i,j,k) = Real(0.5) * (Huon(i,j,k) + Huon(i,j,k-1));
105 FX(i,j,k) = XF(i,j,k) * Real(0.5) * (tke(i-1,j,k,nstp) + tke(i,j,k,nstp) -
106 cff * (grad_ip1 - grad_im1));
107 FXL(i,j,k) = XF(i,j,k) * Real(0.5) * (gls(i-1,j,k,nstp) + gls(i,j,k,nstp) -
108 cff * (gradL_ip1 - gradL_im1));
109 });
110
111 // need EF/FE/FEL from [xlo to xhi ] by [ylo to yhi+1]
112 ParallelFor(grow(ybx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
113 {
114 Real grad_jm1 = (tke(i,j-1,k,nstp) - tke(i,j-2,k,nstp)) * mskv(i,j-1,0);
115 Real grad_jp1 = (tke(i,j+1,k,nstp) - tke(i,j ,k,nstp)) * mskv(i,j+1,0);
116
117 Real gradL_jm1 = (gls(i,j-1,k,nstp) - gls(i,j-2,k,nstp)) * mskv(i,j-1,0);
118 Real gradL_jp1 = (gls(i,j+1,k,nstp) - gls(i,j ,k,nstp)) * mskv(i,j+1,0);
119
120 // Adjust boundaries
121 // TODO: Make sure indices match with what ROMS does
122 if (j == dlo.y-1 && !is_periodic_in_y) {
123 grad_jm1 = tke(i,j,k,nstp) - tke(i,j-1,k,nstp);
124 gradL_jm1 = gls(i,j,k,nstp) - gls(i,j-1,k,nstp);
125 }
126 else if (j == dhi.y+1 && !is_periodic_in_y) {
127 grad_jp1 = tke(i,j,k,nstp) - tke(i,j-1,k,nstp);
128 gradL_jp1 = gls(i,j,k,nstp) - gls(i,j-1,k,nstp);
129 }
130 Real cff = one/Real(6.0);
131 EF(i,j,k) = Real(0.5) * (Hvom(i,j,k) + Hvom(i,j,k-1));
132 FE(i,j,k) = EF(i,j,k) * Real(0.5) * (tke(i,j-1,k,nstp) + tke(i,j,k,nstp) -
133 cff * (grad_jp1 - grad_jm1));
134 FEL(i,j,k) = EF(i,j,k) * Real(0.5) * (gls(i,j-1,k,nstp) + gls(i,j,k,nstp) -
135 cff * (gradL_jp1 - gradL_jm1));
136 });
137
138 Real gamma = one / Real(6.0);
139 Real cff1, cff2, cff3;
140 int indx;
141 // Time step horizontal advection
142 if (iic == ntfirst) {
143 cff1 = one;
144 cff2 = zero;
145 cff3 = Real(0.5) * dt_lev;
146 indx = nstp;
147 } else {
148 cff1 = Real(0.5) + gamma;
149 cff2 = Real(0.5) - gamma;
150 cff3 = (one - gamma) * dt_lev;
151 indx = 1 - nstp;
152 }
153
154 // update tke, gls from [xlo to xhi ] by [ylo to yhi ]
155 // need XF/FX/FXL from [xlo to xhi+1] by [ylo to yhi ]
156 // need EF/FE/FEL from [xlo to xhi ] by [ylo to yhi+1]
157 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
158 {
159 Real cff = Real(0.5) * (Hz(i,j,k) + Hz(i,j,k-1));
160 Real cff4 = cff3 * pm(i,j,0) * pn(i,j,0);
161 Hz_half(i,j,k) = cff - cff4 * (XF(i+1,j,k)-XF(i,j,k)+EF(i,j+1,k)-EF(i,j,k));
162 tke(i,j,k,2) = cff * (cff1*tke(i,j,k,nstp) + cff2*tke(i,j,k,indx)) -
163 cff4 * (FX(i+1,j,k)-FX(i,j,k)+FE(i,j+1,k)-FE(i,j,k));
164 gls(i,j,k,2) = cff * (cff1 * gls(i,j,k,nstp) + cff2 * gls(i,j,k,indx)) -
165 cff4 * (FXL(i+1,j,k)-FXL(i,j,k)+FEL(i,j+1,k)-FEL(i,j,k));
166 tke(i,j,k,nnew) = cff * tke(i,j,k,nstp);
167 gls(i,j,k,nnew) = cff * gls(i,j,k,nstp);
168 });
169
170 // Will do a FillPatch after this, so don't need to do any ghost zones in x,y
171 // Compute vertical advection
172 ParallelFor(convert(bx,IntVect(0,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
173 {
174 // CF and FC/FCL are on rho points
175 CF(i,j,k) = Real(0.5) * (W(i,j,k+1) + W(i,j,k));
176 if (k == 0) {
177 Real cff1_vadv = one / Real(3.0);
178 Real cff2_vadv = Real(5.0) / Real(6.0);
179 Real cff3_vadv = one / Real(6.0);
180 FC(i,j,k) = CF(i,j,k) * (cff1_vadv * tke(i,j,0,nstp) +
181 cff2_vadv * tke(i,j,1,nstp) -
182 cff3_vadv * tke(i,j,2,nstp));
183 FCL(i,j,k) = CF(i,j,k) * (cff1_vadv * gls(i,j,0,nstp) +
184 cff2_vadv * gls(i,j,1,nstp) -
185 cff3_vadv * gls(i,j,2,nstp));
186 } else if (k == N) {
187 Real cff1_vadv = one / Real(3.0);
188 Real cff2_vadv = Real(5.0) / Real(6.0);
189 Real cff3_vadv = one / Real(6.0);
190 FC(i,j,k) = CF(i,j,k) * (cff1_vadv * tke(i,j,k+1, nstp) +
191 cff2_vadv * tke(i,j,k ,nstp)-
192 cff3_vadv * tke(i,j,k-1,nstp));
193 FCL(i,j,k) = CF(i,j,k) * (cff1_vadv * gls(i,j,k+1,nstp) +
194 cff2_vadv * gls(i,j,k ,nstp)-
195 cff3_vadv * gls(i,j,k-1,nstp));
196 } else {
197 Real cff1_vadv = Real(7.0) / Real(12.0);
198 Real cff2_vadv = one / Real(12.0);
199 FC(i,j,k) = CF(i,j,k) * (cff1_vadv * (tke(i,j,k ,nstp) + tke(i,j,k+1,nstp)) -
200 cff2_vadv * (tke(i,j,k-1,nstp) + tke(i,j,k+2,nstp)));
201 FCL(i,j,k) = CF(i,j,k) * (cff1_vadv * (gls(i,j,k ,nstp) + gls(i,j,k+1,nstp)) -
202 cff2_vadv * (gls(i,j,k-1,nstp) + gls(i,j,k+2,nstp)));
203 }
204 });
205
206 // Time-step vertical advection
207 if (iic == ntfirst) {
208 cff3 = Real(0.5) * dt_lev;
209 } else {
210 cff3 = (one - gamma) * dt_lev;
211 }
212 // DO k=1,N-1
213 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
214 {
215 Real cff4 = cff3 * pm(i,j,0) * pn(i,j,0);
216 Hz_half(i,j,k) = Hz_half(i,j,k) - cff4 * (CF(i,j,k)-CF(i,j,k-1));
217 Real cff1_loc = one / Hz_half(i,j,k);
218 tke(i,j,k,2) = cff1_loc * (tke(i,j,k,2) - cff4 * (FC (i,j,k) - FC (i,j,k-1)));
219 gls(i,j,k,2) = cff1_loc * (gls(i,j,k,2) - cff4 * (FCL(i,j,k) - FCL(i,j,k-1)));
220 });
221 }
222
223 for (int icomp=0; icomp<3; icomp++) {
224 FillPatch(lev, t_old[lev], *vec_tke[lev], GetVecOfPtrs(vec_tke), tke_bc(), BdyVars::null, icomp, false, false);
225 FillPatch(lev, t_old[lev], *vec_gls[lev], GetVecOfPtrs(vec_gls), tke_bc(), BdyVars::null, icomp, false, false);
226 }
227}
228
229/**
230 * @param[in ] lev level to operate on
231 * @param[inout] mf_gls turbulent generic length scale
232 * @param[inout] mf_tke turbulent kinetic energy
233 * @param[in ] mf_W vertical velocity
234 * @param[inout] mf_Akv vertical viscosity coefficient
235 * @param[inout] mf_Akt vertical diffusivity coefficients
236 * @param[inout] mf_Akk turbulent kinetic energy vertical diffusion coefficient
237 * @param[inout] mf_Akp turbulent length scale vertical diffusion coefficient
238 * @param[in ] mf_mskr land-sea mask on rho points
239 * @param[in ] mf_msku land-sea mask on u points
240 * @param[in ] mf_mskv land-sea mask on v points
241 * @param[in ] nstp index of last time step in gls and tke MultiFabs
242 * @param[in ] nnew index of time step to update in gls and tke MultiFabs
243 * @param[in ] N number of vertical levels
244 * @param[in ] dt_lev time step at this level
245 */
246void
247REMORA::gls_corrector (int lev, MultiFab* mf_gls, MultiFab* mf_tke,
248 MultiFab& mf_W, MultiFab* mf_Akv, MultiFab* mf_Akt,
249 MultiFab* mf_Akk, MultiFab* mf_Akp,
250 MultiFab* mf_mskr,
251 MultiFab* mf_msku, MultiFab* mf_mskv,
252 const int nstp, const int nnew,
253 const int N, const Real dt_lev)
254{
255 BL_PROFILE("REMORA::gls_corrector()");
256//-----------------------------------------------------------------------
257// Compute several constants.
258//-----------------------------------------------------------------------
259 bool Lmy25 = ((solverChoice.gls_p == zero) &&
260 (solverChoice.gls_n == one) &&
261 (solverChoice.gls_m == one)) ? true : false;
262
263 Real L_sft = vonKar;
264 Real gls_sigp_cb = solverChoice.gls_sigp;
265 Real ogls_sigp = one/gls_sigp_cb;
266
267 Real gls_c3m = solverChoice.gls_c3m;
268 Real gls_c3p = solverChoice.gls_c3p;
269 Real gls_cmu0 = solverChoice.gls_cmu0;
270
271 Real gls_m = solverChoice.gls_m;
272 Real gls_n = solverChoice.gls_n;
273 Real gls_p = solverChoice.gls_p;
274
275 Real gls_Gh0 = solverChoice.gls_Gh0;
276 Real gls_Ghcri = solverChoice.gls_Ghcri;
277 Real gls_Ghmin = solverChoice.gls_Ghmin;
278
279 Real Akv_bak = solverChoice.Akv_bak;
280 Real Akp_bak = solverChoice.Akp_bak;
281 Real Akk_bak = solverChoice.Akk_bak;
282
283 // Akt_bak has one entry per active tracer, so it has to reach the device as an array.
284 // The stratification terms below are specifically the temperature ones, and use
285 // Akt_bak[Temp_comp] to match the Akt component they are paired with.
286 //
287 // A GpuArray captured by value rather than a device allocation: NAT is a compile-time
288 // constant, and a Gpu::DeviceVector here would be freed at function exit while the
289 // kernels launched below could still be reading it -- arena frees are not stream
290 // ordered.
291 GpuArray<Real, NAT> Akt_bak{};
292 for (int n = 0; n < NAT; ++n) {
293 Akt_bak[n] = solverChoice.Akt_bak[n];
294 }
295
296 Real gls_c1 = solverChoice.gls_c1;
297 Real gls_c2 = solverChoice.gls_c2;
298 Real gls_E2 = solverChoice.gls_E2;
299 Real gls_sigk = solverChoice.gls_sigk;
300 auto gls_stability_type = solverChoice.gls_stability_type;
301
302 Real sqrt2 = std::sqrt(two);
303 Real cmu_fac1 = std::pow(solverChoice.gls_cmu0,(-solverChoice.gls_p/solverChoice.gls_n));
304 Real cmu_fac2 = std::pow(solverChoice.gls_cmu0,(Real(3.0)+solverChoice.gls_p/solverChoice.gls_n));
305 Real cmu_fac3 = one/std::pow(solverChoice.gls_cmu0,two);
306
307 Real gls_fac2 = std::pow(solverChoice.gls_cmu0,solverChoice.gls_p)*solverChoice.gls_n*std::pow(vonKar,solverChoice.gls_n);
308 Real gls_fac3 = std::pow(solverChoice.gls_cmu0,solverChoice.gls_p)*solverChoice.gls_n;
309 Real gls_fac4 = std::pow(solverChoice.gls_cmu0,solverChoice.gls_p);
310 Real gls_fac5 = std::pow(Real(0.56),Real(0.5)*solverChoice.gls_n)*std::pow(solverChoice.gls_cmu0,solverChoice.gls_p);
311 Real gls_fac6 = Real(8.0)/std::pow(solverChoice.gls_cmu0,Real(6.0));
312 const bool l_cap = solverChoice.gls_length_cap;
313
314 Real gls_exp1 = one/solverChoice.gls_n;
315 Real tke_exp1 = solverChoice.gls_m/solverChoice.gls_n;
316 Real tke_exp2 = Real(0.5)+solverChoice.gls_m/solverChoice.gls_n;
317 Real tke_exp4 = solverChoice.gls_m+Real(0.5)*solverChoice.gls_n;
318
319 Real cmu0_exp_p = std::pow(gls_cmu0, gls_p);
320 Real gls_cmu0_cube = gls_cmu0 * gls_cmu0 * gls_cmu0;
321
322 Real gls_s0, gls_s1, gls_s2, gls_s4, gls_s5, gls_s6;
323 Real gls_b0, gls_b1, gls_b2, gls_b3, gls_b4, gls_b5;
324 Real my_Sm2, my_Sh1, my_Sh2, my_Sm3, my_Sm4;
325
326 // Compute parameters for Canuto et al. (2001) stability functions.
327 // (Canuto, V.M., Cheng, H.Y., and Dubovikov, M.S., 2001: Ocean
328 // turbulence. Part I: One-point closure model - momentum and
329 // heat vertical diffusivities, JPO, 1413-1426).
330
331 if (solverChoice.gls_stability_type == GLS_StabilityType::Canuto_A ||
332 solverChoice.gls_stability_type == GLS_StabilityType::Canuto_B) {
333
334 gls_s0=Real(3.0)/Real(2.0)*solverChoice.gls_L1*solverChoice.gls_L5*solverChoice.gls_L5;
335 gls_s1=-solverChoice.gls_L4*(solverChoice.gls_L6+solverChoice.gls_L7)
336 +two*solverChoice.gls_L4*solverChoice.gls_L5*
337 (solverChoice.gls_L1-one/Real(3.0)*solverChoice.gls_L2-solverChoice.gls_L3)
338 +Real(3.0)/Real(2.0)*
339 solverChoice.gls_L1*solverChoice.gls_L5*solverChoice.gls_L8;
340 gls_s2=Real(-3.0)/Real(8.0)*solverChoice.gls_L1
341 *(solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
342 gls_s4=two*solverChoice.gls_L5;
343 gls_s5=two*solverChoice.gls_L4;
344 gls_s6=two/Real(3.0)*solverChoice.gls_L5
345 *(Real(3.0)*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2)-
346 Real(0.5)*solverChoice.gls_L5*solverChoice.gls_L1*(Real(3.0)*solverChoice.gls_L3-solverChoice.gls_L2)+
347 Real(3.0)/Real(4.0)*solverChoice.gls_L1*(solverChoice.gls_L6-solverChoice.gls_L7);
348 gls_b0=Real(3.0)*solverChoice.gls_L5*solverChoice.gls_L5;
349 gls_b1=solverChoice.gls_L5*(Real(7.0)*solverChoice.gls_L4+Real(3.0)*solverChoice.gls_L8);
350 gls_b2=solverChoice.gls_L5*solverChoice.gls_L5*(Real(3.0)*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2)-
351 Real(3.0)/Real(4.0)*(solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
352 gls_b3=solverChoice.gls_L4*(Real(4.0)*solverChoice.gls_L4+Real(3.0)*solverChoice.gls_L8);
353 gls_b5=one/Real(4.0)*(solverChoice.gls_L2*solverChoice.gls_L2-Real(3.0)*solverChoice.gls_L3*solverChoice.gls_L3)*
354 (solverChoice.gls_L6*solverChoice.gls_L6-solverChoice.gls_L7*solverChoice.gls_L7);
355 gls_b4=solverChoice.gls_L4*(solverChoice.gls_L2*solverChoice.gls_L6-Real(3.0)*solverChoice.gls_L3*solverChoice.gls_L7-
356 solverChoice.gls_L5*(solverChoice.gls_L2*solverChoice.gls_L2-solverChoice.gls_L3*solverChoice.gls_L3))+solverChoice.gls_L5*solverChoice.gls_L8*
357 (Real(3.0)*solverChoice.gls_L3*solverChoice.gls_L3-solverChoice.gls_L2*solverChoice.gls_L2);
358 my_Sm2 = zero;
359 my_Sm3 = zero;
360 my_Sm4 = zero;
361 my_Sh1 = zero;
362 my_Sh2 = zero;
363 } else {
364 gls_s0 = zero;
365 gls_s1 = zero;
366 gls_s2 = zero;
367 gls_s4 = zero;
368 gls_s5 = zero;
369 gls_s6 = zero;
370 gls_b0 = zero;
371 gls_b1 = zero;
372 gls_b2 = zero;
373 gls_b3 = zero;
374 gls_b4 = zero;
375 gls_b5 = zero;
376 my_Sm2=Real(9.0)*solverChoice.my_A1*solverChoice.my_A2;
377 my_Sm3=solverChoice.my_A1*(one-Real(3.0)*solverChoice.my_C1-Real(6.0)*solverChoice.my_A1/solverChoice.my_B1);
378 my_Sm4=Real(18.0)*solverChoice.my_A1*solverChoice.my_A1+Real(9.0)*solverChoice.my_A1*solverChoice.my_A2;
379 my_Sh1=solverChoice.my_A2*(one-Real(6.0)*solverChoice.my_A1/solverChoice.my_B1);
380 my_Sh2=Real(3.0)*solverChoice.my_A2*(Real(6.0)*solverChoice.my_A1+solverChoice.my_B2);
381 }
382
383 Real Zos_min = std::max(solverChoice.Zos, Real(0.0001));
384 Real Zos_eff = Zos_min;
385 Real Gadv = one/Real(3.0);
386 Real eps = Real(1.0e-10);
387
388 const BoxArray& ba = cons_old[lev]->boxArray();
389 const DistributionMapping& dm = cons_old[lev]->DistributionMap();
390
391 int ncomp_w = 0;
392 int dU_comp = ncomp_w++;
393 int dV_comp = ncomp_w++;
394 int CF_comp = ncomp_w++;
395
396 int ncomp = 0;
397 int shear2_comp = ncomp++;
398 int shear2_cache_comp = ncomp++;
399 int buoy2_comp = ncomp++;
400
401 MultiFab mf_w(convert(ba, IntVect(0,0,1)),dm,ncomp_w,IntVect(NGROW,NGROW,0));
402 MultiFab mf(ba,dm,ncomp,IntVect(NGROW,NGROW,0));
403
404 const Box& domain = geom[0].Domain();
405 const auto dlo = amrex::lbound(domain);
406 const auto dhi = amrex::ubound(domain);
407
408 GeometryData const& geomdata = geom[0].data();
409 bool is_periodic_in_x = geomdata.isPeriodic(0);
410 bool is_periodic_in_y = geomdata.isPeriodic(1);
411
412 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi )
413 {
414 Box bx = mfi.tilebox();
415 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
416
417 Box bxD = bx;
418 bxD.makeSlab(2,0);
419 Box gbx1D = gbx1;
420 gbx1D.makeSlab(2,0);
421
422 Array4<Real> const& Hz = vec_Hz[lev]->array(mfi);
423 Array4<Real> const& u = xvel_old[lev]->array(mfi);
424 Array4<Real> const& v = yvel_old[lev]->array(mfi);
425
426 auto dU = mf_w.array(mfi,dU_comp);
427 auto dV = mf_w.array(mfi,dV_comp);
428 auto CF = mf_w.array(mfi,CF_comp);
429 auto shear2_cached = mf.array(mfi,shear2_cache_comp);
430
431 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
432 {
433 CF(i,j,0) = zero;
434 dU(i,j,0) = zero;
435 dV(i,j,0) = zero;
436 for (int k=1; k<=N; k++) {
437 Real cff = one / (two * Hz(i,j,k) + Hz(i,j,k-1)*(two - CF(i,j,k-1)));
438 CF(i,j,k) = cff * Hz(i,j,k);
439 dU(i,j,k)=cff*(Real(3.0)*(u(i ,j,k)-u(i, j,k-1)+
440 u(i+1,j,k)-u(i+1,j,k-1))-Hz(i,j,k-1)*dU(i,j,k-1));
441 dV(i,j,k)=cff*(Real(3.0)*(v(i,j ,k)-v(i,j ,k-1)+
442 v(i,j+1,k)-v(i,j+1,k-1))-Hz(i,j,k-1)*dV(i,j,k-1));
443 }
444 dU(i,j,N+1) = zero;
445 dV(i,j,N+1) = zero;
446 for (int k=N; k>=1; k--) {
447 dU(i,j,k) = dU(i,j,k) - CF(i,j,k) * dU(i,j,k+1);
448 dV(i,j,k) = dV(i,j,k) - CF(i,j,k) * dV(i,j,k+1);
449 }
450 shear2_cached(i,j,0) = zero;
451 for (int k=1; k<=N; k++) {
452 shear2_cached(i,j,k) = dU(i,j,k) * dU(i,j,k) + dV(i,j,k) * dV(i,j,k);
453 }
454 });
455 }
456
457 // While potentially counterintuitive, this is what ROMS does for handling shear2 at all boundaries, even
458 // periodic
459 (*physbcs[lev])(mf,*mf_mskr,shear2_cache_comp,1,mf.nGrowVect(),t_new[lev],foextrap_bc());
460 mf.setVal(zero,CF_comp,1);
461
462 int ncomp_fab = 0;
463 int tmp_buoy_comp = ncomp_fab++;
464 int tmp_shear_comp = ncomp_fab++;
465 int curvK_comp = ncomp_fab++;
466 int curvP_comp = ncomp_fab++;
467 int FXK_comp = ncomp_fab++;
468 int FXP_comp = ncomp_fab++;
469 int FEK_comp = ncomp_fab++;
470 int FEP_comp = ncomp_fab++;
471 int FCK_comp = ncomp_fab++;
472 int FCP_comp = ncomp_fab++;
473 int BCK_comp = ncomp_fab++;
474 int BCP_comp = ncomp_fab++;
475
476 for ( MFIter mfi(*mf_gls, TilingIfNotGPU()); mfi.isValid(); ++mfi )
477 {
478 Box bx = mfi.tilebox();
479 Box xbx = surroundingNodes(bx,0);
480 Box ybx = surroundingNodes(bx,1);
481 Box gbx1 = grow(bx,IntVect(NGROW-1,NGROW-1,0));
482
483 Box bx_rho = bx;
484 bx_rho.convert(IntVect(0,0,0));
485 Box bx_growloxy = growLo(growLo(grow(bx,IntVect(0,0,-1)),0,1),1,1);
486
487 Box bxD = bx;
488 bxD.makeSlab(2,0);
489 Box gbx1D = gbx1;
490 gbx1D.makeSlab(2,0);
491
492 int ncompbc = 1;
493 Vector<BCRec> bcrs_x(ncompbc);
494 Vector<BCRec> bcrs_y(ncompbc);
495 amrex::setBC(xbx,domain,xvel_bc(),0,1,domain_bcs_type,bcrs_x);
496 amrex::setBC(ybx,domain,yvel_bc(),0,1,domain_bcs_type,bcrs_y);
497
498 Array4<Real const> const& W = mf_W.const_array(mfi);
499 Array4<Real> const& Hz = vec_Hz[lev]->array(mfi);
500 Array4<Real> const& pm = vec_pm[lev]->array(mfi);
501 Array4<Real> const& pn = vec_pn[lev]->array(mfi);
502 Array4<Real> const& Lscale = vec_Lscale[lev]->array(mfi);
503
504 Array4<Real> const& Huon = vec_Huon[lev]->array(mfi);
505 Array4<Real> const& Hvom = vec_Hvom[lev]->array(mfi);
506 Array4<Real> const& z_w = vec_z_w[lev]->array(mfi);
507
508 Array4<Real> const& tke = mf_tke->array(mfi);
509 Array4<Real> const& gls = mf_gls->array(mfi);
510
511 Array4<Real const> const& sustr = vec_sustr[lev]->const_array(mfi);
512 Array4<Real const> const& svstr = vec_svstr[lev]->const_array(mfi);
513 Array4<Real const> const& bustr = vec_bustr[lev]->const_array(mfi);
514 Array4<Real const> const& bvstr = vec_bvstr[lev]->const_array(mfi);
515 Array4<Real const> const& msku = mf_msku->const_array(mfi);
516 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
517
518 Array4<Real> const& ZoBot = vec_ZoBot[lev]->array(mfi);
519
520 FArrayBox fab(gbx1,ncomp_fab, amrex::The_Async_Arena()); fab.template setVal<RunOn::Device>(zero);
521
522 auto CF = mf_w.array(mfi,CF_comp);
523 auto shear2 = mf.array(mfi,shear2_comp);
524 auto shear2_cached = mf.array(mfi,shear2_cache_comp);
525 auto buoy2 = mf.array(mfi,buoy2_comp);
526 Array4<Real> const& bvf = vec_bvf[lev]->array(mfi);
527
528 auto tmp_buoy = fab.array(tmp_buoy_comp);
529 auto tmp_shear = fab.array(tmp_shear_comp);
530 auto curvK = fab.array(curvK_comp);
531 auto curvP = fab.array(curvP_comp);
532 auto FXK = fab.array(FXK_comp);
533 auto FXP = fab.array(FXP_comp);
534 auto FEK = fab.array(FEK_comp);
535 auto FEP = fab.array(FEP_comp);
536 auto FCK = fab.array(FCK_comp);
537 auto FCP = fab.array(FCP_comp);
538 auto BCK = fab.array(BCK_comp);
539 auto BCP = fab.array(BCP_comp);
540
541 auto Akt = mf_Akt->array(mfi);
542 auto Akv = mf_Akv->array(mfi);
543 auto Akp = mf_Akp->array(mfi);
544 auto Akk = mf_Akk->array(mfi);
545
546 ParallelFor(bx_growloxy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
547 {
548 tmp_buoy(i,j,k)=Real(0.25) * (bvf(i,j,k) + bvf(i+1,j,k) + bvf(i,j+1,k)+bvf(i+1,j+1,k));
549 tmp_shear(i,j,k)=Real(0.25) * (shear2_cached(i,j,k) + shear2_cached(i+1,j,k) + shear2_cached(i,j+1,k)+shear2_cached(i+1,j+1,k));
550 });
551
552 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
553 {
554 buoy2(i,j,k)=Real(0.25) * (tmp_buoy(i,j,k) + tmp_buoy(i-1,j,k) + tmp_buoy(i,j-1,k)+tmp_buoy(i-1,j-1,k));
555 shear2(i,j,k)=Real(0.25) * (tmp_shear(i,j,k) + tmp_shear(i-1,j,k) + tmp_shear(i,j-1,k)+tmp_shear(i-1,j-1,k));
556 });
557
558 //Time step advective terms
559 ParallelFor(growLo(grow(xbx,IntVect(0,0,-1)),0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
560 {
561 Real gradK, gradK_ip1, gradP, gradP_ip1;
562
563 if (i == dlo.x-1 && !is_periodic_in_x) {
564 gradK_ip1 = tke(i+1,j,k,2)-tke(i ,j,k,2);
565 gradK = gradK_ip1;
566 gradP_ip1 = gls(i+1,j,k,2)-gls(i ,j,k,2);
567 gradP = gradP_ip1;
568 } else if (i == dhi.x+1 && !is_periodic_in_x) {
569 gradK = tke(i ,j,k,2)-tke(i-1,j,k,2);
570 gradK_ip1 = gradK;
571 gradP = gls(i ,j,k,2)-gls(i-1,j,k,2);
572 gradP_ip1 = gradP;
573 } else {
574 gradK = (tke(i ,j,k,2)-tke(i-1,j,k,2)) * msku(i ,j,0);
575 gradK_ip1 = (tke(i+1,j,k,2)-tke(i ,j,k,2)) * msku(i+1,j,0);
576 gradP = (gls(i ,j,k,2)-gls(i-1,j,k,2)) * msku(i ,j,0);
577 gradP_ip1 = (gls(i+1,j,k,2)-gls(i ,j,k,2)) * msku(i+1,j,0);
578 }
579
580 curvK(i,j,k) = gradK_ip1 - gradK;
581 curvP(i,j,k) = gradP_ip1 - gradP;
582 });
583 ParallelFor(grow(xbx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
584 {
585 Real cff = Real(0.5) * (Huon(i,j,k) + Huon(i,j,k-1));
586 Real cff1 = (cff > zero) ? curvK(i-1,j,k) : curvK(i,j,k);
587 Real cff2 = (cff > zero) ? curvP(i-1,j,k) : curvP(i,j,k);
588
589 FXK(i,j,k) = cff * Real(0.5) * (tke(i-1,j,k,2)+tke(i,j,k,2)-Gadv*cff1);
590 FXP(i,j,k) = cff * Real(0.5) * (gls(i-1,j,k,2)+gls(i,j,k,2)-Gadv*cff2);
591 });
592
593 //Time step advective terms
594 ParallelFor(growLo(grow(ybx,IntVect(0,0,-1)),1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
595 {
596 Real gradK = (tke(i,j ,k,2)-tke(i,j-1,k,2)) * mskv(i,j ,0);
597 Real gradK_jp1 = (tke(i,j+1,k,2)-tke(i,j ,k,2)) * mskv(i,j+1,0);
598 Real gradP = (gls(i,j ,k,2)-gls(i,j-1,k,2)) * mskv(i,j ,0);
599 Real gradP_jp1 = (gls(i,j+1,k,2)-gls(i,j ,k,2)) * mskv(i,j+1,0);
600
601 if (j == dlo.y-1 && !is_periodic_in_y) {
602 gradK = gradK_jp1;
603 gradP = gradP_jp1;
604 }
605 else if (j == dhi.y+1 && !is_periodic_in_y) {
606 gradK_jp1 = gradK;
607 gradP_jp1 = gradP;
608 }
609
610 curvK(i,j,k) = gradK_jp1 - gradK;
611 curvP(i,j,k) = gradP_jp1 - gradP;
612 });
613 ParallelFor(grow(ybx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
614 {
615 Real cff = Real(0.5) * (Hvom(i,j,k) + Hvom(i,j,k-1));
616 Real cff1 = (cff > zero) ? curvK(i,j-1,k) : curvK(i,j,k);
617 Real cff2 = (cff > zero) ? curvP(i,j-1,k) : curvP(i,j,k);
618
619 FEK(i,j,k) = cff * Real(0.5) * (tke(i,j-1,k,2)+tke(i,j,k,2)-Gadv*cff1);
620 FEP(i,j,k) = cff * Real(0.5) * (gls(i,j-1,k,2)+gls(i,j,k,2)-Gadv*cff2);
621 });
622
623 Real gls_Kmin = solverChoice.gls_Kmin;
624 Real gls_Pmin = solverChoice.gls_Pmin;
625 ParallelFor(grow(bx,IntVect(0,0,-1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
626 {
627 Real cff = dt_lev * pm(i,j,0) * pn(i,j,0);
628 tke(i,j,k,nnew) = tke(i,j,k,nnew) - cff * (FXK(i+1,j ,k)-FXK(i,j,k)+
629 FEK(i ,j+1,k)-FEK(i,j,k));
630 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew), gls_Kmin);
631
632 gls(i,j,k,nnew) = gls(i,j,k,nnew) - cff * (FXP(i+1,j ,k)-FXP(i,j,k)+
633 FEP(i ,j+1,k)-FEP(i,j,k));
634 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew), gls_Pmin);
635 });
636
637 // Vertical advection
638 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
639 {
640 Real cff1 = Real(7.0) / Real(12.0);
641 Real cff2 = one / Real(12.0);
642 for (int k=1; k<=N-1; k++) {
643 Real cff = Real(0.5) * (W(i,j,k+1)+W(i,j,k));
644 FCK(i,j,k) = cff * (cff1 * (tke(i,j,k ,2)+tke(i,j,k+1,2))-
645 cff2 * (tke(i,j,k-1,2)+tke(i,j,k+2,2)));
646 FCP(i,j,k) = cff * (cff1 * (gls(i,j,k ,2)+gls(i,j,k+1,2))-
647 cff2 * (gls(i,j,k-1,2)+gls(i,j,k+2,2)));
648 }
649 cff1 = one/Real(3.0);
650 cff2 = Real(5.0)/Real(6.0);
651 Real cff3 = one / Real(6.0);
652 Real cff = Real(0.5) * (W(i,j,0)+W(i,j,1));
653 FCK(i,j,0) = cff * (cff1 * tke(i,j,0,2)+cff2 * tke(i,j,1,2)-cff3 * tke(i,j,2,2));
654 FCP(i,j,0) = cff * (cff1 * gls(i,j,0,2)+cff2 * gls(i,j,1,2)-cff3 * gls(i,j,2,2));
655 cff = Real(0.5) * (W(i,j,N+1)+W(i,j,N));
656 FCK(i,j,N) = cff * (cff1 * tke(i,j,N+1,2)+cff2*tke(i,j,N,2)-cff3*tke(i,j,N-1,2));
657 FCP(i,j,N) = cff * (cff1 * gls(i,j,N+1,2)+cff2*gls(i,j,N,2)-cff3*gls(i,j,N-1,2));
658 });
659 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
660 {
661 Real cff = dt_lev * pm(i,j,0) * pn(i,j,0);
662 tke(i,j,k,nnew) = tke(i,j,k,nnew) - cff*(FCK(i,j,k )-FCK(i,j,k-1));
663 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew),gls_Kmin);
664 gls(i,j,k,nnew) = gls(i,j,k,nnew) - cff*(FCP(i,j,k )-FCP(i,j,k-1));
665 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew),gls_Pmin);
666 });
667
668 // Compute vertical mixing, turbulent production and turbulent
669 // dissipation.
670 //
671 Real cff = -Real(0.5) * dt_lev;
672 ParallelFor(convert(bx,IntVect(0,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
673 {
674 if (k==0 or k==N) {
675 FCK(i,j,k) = zero;
676 FCP(i,j,k) = zero;
677 } else {
678 FCK(i,j,k) = cff * (Akk(i,j,k) + Akk(i,j,k+1)) / Hz(i,j,k);
679 FCP(i,j,k) = cff * (Akp(i,j,k) + Akp(i,j,k+1)) / Hz(i,j,k);
680 }
681 });
682 // Compute production and dissipation terms.
683 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
684 {
685 // Compute shear and buoyant production of turbulent energy (m3/s3)
686 // at W-points (ignore small negative values of buoyancy).
687 Real strat2 = buoy2(i,j,k);
688 Real gls_c3 = (strat2 > zero) ? gls_c3m : gls_c3p;
689 Real Kprod = shear2(i,j,k) * (Akv(i,j,k)-Akv_bak) -
690 strat2 * (Akt(i,j,k,Temp_comp)-Akt_bak[Temp_comp]);
691 Real Pprod = gls_c1 * shear2(i,j,k) * (Akv(i,j,k)-Akv_bak) -
692 gls_c3 * strat2 * (Akt(i,j,k,Temp_comp)-Akt_bak[Temp_comp]);
693
694 // If negative production terms, then add buoyancy to dissipation terms
695 // (BCK and BCP) below, using "cff1" and "cff2" as the on/off switch.
696 Real cff1 = (Kprod < zero) ? zero : one;
697 Real cff2 = (Pprod < zero) ? zero : one;
698 Kprod = (Kprod < zero) ? Kprod + strat2*(Akt(i,j,k,Temp_comp)-Akt_bak[Temp_comp]) : Kprod;
699 Pprod = (Pprod < zero) ? Pprod + gls_c3*strat2*(Akt(i,j,k,Temp_comp)-Akt_bak[Temp_comp]) : Pprod;
700 // Time-step shear and buoyancy production terms.
701 Real cff_Hz = Real(0.5) * (Hz(i,j,k) + Hz(i,j,k-1));
702 tke(i,j,k,nnew) = tke(i,j,k,nnew)+dt_lev * cff_Hz * Kprod;
703 gls(i,j,k,nnew) = gls(i,j,k,nnew)+dt_lev
704 *cff_Hz*Pprod*gls(i,j,k,nstp) / std::max(tke(i,j,k,nstp),gls_Kmin);
705
706 Real gls_exp_exp1 = std::pow(gls(i,j,k,nstp),gls_exp1);
707 Real gls_exp_mexp1 = one / (gls_exp_exp1);
708 Real tke_exp_mexp1 = std::pow(tke(i,j,k,nstp),-tke_exp1);
709 Real tke_exp_exp2 = std::pow(tke(i,j,k,nstp),tke_exp2);
710
711 // Compute dissipation of turbulent energy (m3/s3).
712 Real wall_fac = one;
713 if (Lmy25) {
714 wall_fac=one+gls_E2/(vonKar*vonKar)*
715 std::pow(gls_exp_exp1*cmu_fac1*
716 tke_exp_mexp1*
717 (one/ (z_w(i,j,k)-z_w(i,j,0))),2)+
718 Real(0.25)/(vonKar*vonKar)*
719 std::pow(gls_exp_exp1*cmu_fac1*
720 tke_exp_mexp1*
721 (one/ (z_w(i,j,N+1)-z_w(i,j,k))),2);
722 }
723 BCK(i,j,k)=cff_Hz*(one+dt_lev*
724 gls_exp_mexp1*cmu_fac2*
725 tke_exp_exp2+
726 dt_lev*(one-cff1)*strat2*
727 (Akt(i,j,k,Temp_comp)-Akt_bak[Temp_comp])/
728 tke(i,j,k,nstp))-
729 FCK(i,j,k)-FCK(i,j,k-1);
730 BCP(i,j,k)=cff_Hz*(one+dt_lev*gls_c2*wall_fac*
731 gls_exp_mexp1*cmu_fac2*
732 tke_exp_exp2+
733 dt_lev*(one-cff2)*gls_c3*strat2*
734 (Akt(i,j,k,Temp_comp)-Akt_bak[Temp_comp])/
735 tke(i,j,k,nstp))-
736 FCP(i,j,k)-FCP(i,j,k-1);
737 });
738
739 // Compute production and dissipation terms.
740 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
741 {
742 Real Zob_min = std::max(ZoBot(i,j,0), Real(0.0001));
743 //----------------------------------------------------------------------
744 // Time-step dissipation and vertical diffusion terms implicitly.
745 //----------------------------------------------------------------------
746 //
747 // Set Dirichlet surface and bottom boundary conditions. Compute
748 // surface roughness from wind stress (Charnok) and set Craig and
749 // Banner wave breaking surface flux, if appropriate.
750
751
752 tke(i,j,N+1,nnew)=std::max(cmu_fac3*Real(0.5)*
753 std::sqrt((sustr(i,j,0)+sustr(i+1,j,0))*(sustr(i,j,0)+sustr(i+1,j,0))+
754 (svstr(i,j,0)+svstr(i,j+1,0))*(svstr(i,j,0)+svstr(i,j+1,0))),
755 gls_Kmin);
756 tke(i,j,0,nnew)=std::max(cmu_fac3*Real(0.5)*
757 std::sqrt((bustr(i,j,0)+bustr(i+1,j,0))*(bustr(i,j,0)+bustr(i+1,j,0))+
758 (bvstr(i,j,0)+bvstr(i,j+1,0))*(bvstr(i,j,0)+bvstr(i,j+1,0))),
759 gls_Kmin);
760
761 gls(i,j,N+1,nnew)=std::max(cmu0_exp_p*
762 std::pow(tke(i,j,N+1,nnew),gls_m)*
763 std::pow(L_sft*Zos_eff,gls_n), gls_Pmin);
764 Real cff_gls = gls_fac4*std::pow(vonKar*Zob_min,gls_n);
765 gls(i,j,0,nnew)=std::max(cff_gls*std::pow(tke(i,j,0,nnew),(gls_m)), gls_Pmin);
766
767 // Solve tri-diagonal system for turbulent kinetic energy.
768 // Might be N instead of N-1?
769 Real tke_fluxt = zero;
770 Real tke_fluxb = zero;
771 Real cff_BCK = one/BCK(i,j,N);
772 CF(i,j,N)=cff_BCK*FCK(i,j,N-1);
773 tke(i,j,N,nnew)=cff_BCK*(tke(i,j,N,nnew)+tke_fluxt);
774 for (int k=N-1;k>=1;k--) {
775 cff_BCK = one / (BCK(i,j,k)-CF(i,j,k+1)*FCK(i,j,k));
776 CF(i,j,k) = cff_BCK * FCK(i,j,k-1);
777 tke(i,j,k,nnew) = cff_BCK * (tke(i,j,k,nnew) - FCK(i,j,k) * tke(i,j,k+1,nnew));
778 }
779 tke(i,j,1,nnew) = tke(i,j,1,nnew) - cff_BCK * tke_fluxb;
780 tke(i,j,1,nnew) = std::max(tke(i,j,1,nnew),gls_Kmin);
781 for (int k=2;k<=N;k++) {
782 tke(i,j,k,nnew) = tke(i,j,k,nnew) - CF(i,j,k) * tke(i,j,k-1,nnew);
783 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew), gls_Kmin);
784 }
785
786 // Solve tri-diagonal system for generic statistical field.
787 Real cff_tke = Real(0.5) * (tke(i,j,N+1,nnew) + tke(i,j,N,nnew));
788 Real gls_fluxt = dt_lev*gls_fac3*std::pow(cff_tke,gls_m)*
789 std::pow(L_sft,(gls_n))*
790 std::pow(Zos_eff+Real(0.5)*Hz(i,j,N),gls_n-one)*
791 Real(0.5)*(Akp(i,j,N+1)+Akp(i,j,N));
792 cff_tke=Real(0.5)*(tke(i,j,0,nnew)+tke(i,j,1,nnew));
793 Real gls_fluxb = dt_lev*gls_fac2*std::pow(cff_tke,gls_m)*
794 std::pow(Real(0.5)*Hz(i,j,0)+Zob_min,gls_n-one)*
795 Real(0.5)*(Akp(i,j,0)+Akp(i,j,1));
796 Real cff_BCP = one / BCP(i,j,N);
797 CF(i,j,N) = cff_BCP * FCP(i,j,N-1);
798 gls(i,j,N,nnew)=cff_BCP*(gls(i,j,N,nnew)-gls_fluxt);
799 for (int k=N-1;k>=1;k--) {
800 cff_BCP = one / (BCP(i,j,k)-CF(i,j,k+1)*FCP(i,j,k));
801 CF(i,j,k) = cff_BCP * FCP(i,j,k-1);
802 gls(i,j,k,nnew) = cff_BCP * (gls(i,j,k,nnew) - FCP(i,j,k)*gls(i,j,k+1,nnew));
803 }
804 gls(i,j,1,nnew) = gls(i,j,1,nnew)-cff_BCP*gls_fluxb;
805 for (int k=2; k<=N; k++) {
806 gls(i,j,k,nnew) = gls(i,j,k,nnew) - CF(i,j,k) * gls(i,j,k-1,nnew);
807 }
808 });
809
810 // Compute vertical mixing coefficients (m2/s).
811 ParallelFor(grow(bx,2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
812 {
813 tke(i,j,k,nnew) = std::max(tke(i,j,k,nnew),gls_Kmin);
814 gls(i,j,k,nnew) = std::max(gls(i,j,k,nnew),gls_Pmin);
815 Real gls_comparison = gls_fac5 *
816 std::pow(tke(i,j,k,nnew),tke_exp4)*
817 std::pow(std::sqrt(std::max(Real(0.0),
818 buoy2(i,j,k)))+eps,-gls_n);
819 gls(i,j,k,nnew) = (gls_n >= Real(0.0)) ? std::min(gls(i,j,k,nnew),gls_comparison) : std::max(gls(i,j,k,nnew),gls_comparison);
820 Real Ls_lmt;
821 Real Ls_unlmt=std::max(eps,
822 std::pow(gls(i,j,k,nnew),( gls_exp1))*cmu_fac1*
823 std::pow(tke(i,j,k,nnew),(-tke_exp1)));
824 // Some problems are very sensitive to this condition (ultimate cause of
825 // some discrepancies in BoundaryLayer test between CPU and GPU)
826 Ls_lmt = (buoy2(i,j,k) > Real(0.0)) ? std::min(Ls_unlmt,
827 std::sqrt(Real(0.56)*tke(i,j,k,nnew)/
828 (std::max(Real(0.0),buoy2(i,j,k))+eps))) : Ls_unlmt;
829 //
830 // L <= kappa*Depth (ROMS gls_corstep.F). The Galperin limit above
831 // applies only where buoy2 > 0; this bounds L everywhere else.
832 if (l_cap) { Ls_lmt = std::min(Ls_lmt, vonKar * (z_w(i,j,N+1) - z_w(i,j,0))); }
833 //
834 // Recompute gls based on limited length scale
835 //
836 gls(i,j,k,nnew)=std::max(cmu0_exp_p*
837 std::pow(tke(i,j,k,nnew),gls_m)*
838 std::pow(Ls_lmt,gls_n), gls_Pmin);
839
840 // Compute nondimensional stability functions for tracers (Sh) and
841 // momentum (Sm).
842 Real Sh, Sm;
843 Real Gh=std::min(gls_Gh0,-buoy2(i,j,k)*Ls_lmt*Ls_lmt/
844 (two*tke(i,j,k,nnew)));
845 Gh=std::min(Gh,Gh-(Gh-gls_Ghcri)*(Gh-gls_Ghcri)/
846 (Gh+gls_Gh0-two*gls_Ghcri));
847 Gh=std::max(Gh,gls_Ghmin);
848
849 if (gls_stability_type == GLS_StabilityType::Canuto_A ||
850 gls_stability_type == GLS_StabilityType::Canuto_B) {
851 //
852 // Canuto stability: Compute shear number.
853 //
854 Real Gm=(gls_b0/gls_fac6-gls_b1*Gh+gls_b3*gls_fac6*(Gh*Gh))/
855 (gls_b2-gls_b4*gls_fac6*Gh);
856 Gm=std::min(Gm,shear2(i,j,k)*Ls_lmt*Ls_lmt/
857 (two*tke(i,j,k,nnew)));
858 /////Gm=std::min(Gm,(gls_s1*gls_fac6*Gh-gls_s0)/(gls_s2*gls_fac6));
859 //
860 // Compute stability functions
861 //
862 Real stab_cff=gls_b0-gls_b1*gls_fac6*Gh+gls_b2*gls_fac6*Gm+
863 gls_b3*gls_fac6*gls_fac6*Gh*Gh-gls_b4*gls_fac6*gls_fac6*Gh*Gm+
864 gls_b5*gls_fac6*gls_fac6*Gm*Gm;
865 Sm=(gls_s0-gls_s1*gls_fac6*Gh+gls_s2*gls_fac6*Gm)/stab_cff;
866 Sh=(gls_s4-gls_s5*gls_fac6*Gh+gls_s6*gls_fac6*Gm)/stab_cff;
867 Sm=std::max(Sm,Real(0.0));
868 Sh=std::max(Sh,Real(0.0));
869
870 //
871 // Relate Canuto stability to ROMS notation
872 //
873 Sm=Sm*sqrt2/(gls_cmu0_cube);
874 Sh=Sh*sqrt2/gls_cmu0_cube;
875 } else if (gls_stability_type == GLS_StabilityType::Galperin) {
876 Real cff_galperin = one - my_Sh2*Gh;
877 Sh = my_Sh1 / cff_galperin;
878 Sm = (my_Sm3+Sh*Gh*my_Sm4)/(one-my_Sm2*Gh);
879 } else {
880 Sh = zero;
881 Sm = zero;
882 }
883
884 // Compute vertical mixing (m2/s) coefficients of momentum and
885 // tracers. Average ql over the two timesteps rather than using
886 // the new Lscale and just averaging tke.
887
888 Real ql=sqrt2*Real(0.5)*(Ls_lmt*std::sqrt(tke(i,j,k,nnew))+
889 Lscale(i,j,k)*std::sqrt(tke(i,j,k,nstp)));
890 Akv(i,j,k)=Akv_bak+Sm*ql;
891 for (int n=0; n<NAT; n++) {
892 Akt(i,j,k,n)=Akt_bak[n]+Sh*ql;
893 }
894
895 // Compute vertical mixing (m2/s) coefficients of turbulent kinetic
896 // energy and generic statistical field.
897
898 Akk(i,j,k)=Akk_bak+Sm*ql/gls_sigk;
899 Akp(i,j,k)=Akp_bak+Sm*ql*ogls_sigp;
900
901 // Save limited length scale.
902 Lscale(i,j,k)=Ls_lmt;
903 });
904
905 // Uncomment to print where the cap binds. CPU reduction every step --
906 // do not leave enabled on a GPU build.
907 // {
908 // Real dmax = Real(0.0), dcap = Real(0.0);
909 // amrex::LoopOnCpu(grow(bx,2,-1), [&] (int i, int j, int k) noexcept {
910 // Real Ls_u = std::max(eps,
911 // std::pow(gls(i,j,k,nnew),gls_exp1)*cmu_fac1*
912 // std::pow(tke(i,j,k,nnew),-tke_exp1));
913 // dmax = std::max(dmax, Ls_u);
914 // dcap = std::max(dcap, vonKar*(z_w(i,j,N+1)-z_w(i,j,0)));
915 // });
916 // amrex::Print() << " [LSCALE] max uncapped L = " << dmax
917 // << " m, cap = " << dcap << " m\n";
918 // }
919
920 ParallelFor(bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
921 {
922 Real Zob_min = std::max(ZoBot(i,j,0), Real(0.0001));
923 Akv(i,j,N+1)=Akv_bak+L_sft*Zos_eff*gls_cmu0*
924 std::sqrt(tke(i,j,N+1,nnew));
925 Akv(i,j,0)=Akv_bak+vonKar*Zob_min*gls_cmu0*
926 std::sqrt(tke(i,j,0,nnew));
927
928 Akk(i,j,N+1)=Akk_bak+Akv(i,j,N+1)/gls_sigk;
929 Akk(i,j,0)=Akk_bak+Akv(i,j,0)/gls_sigk;
930 Akp(i,j,N+1)=Akp_bak+Akv(i,j,N+1)*ogls_sigp;
931 Akp(i,j,0)=Akp_bak+Akv(i,j,0)/gls_sigp_cb;
932
933 for (int n=0; n<NAT; n++) {
934 Akt(i,j,N+1,n) = Akt_bak[n];
935 Akt(i,j,0,n) = Akt_bak[n];
936 }
937 });
938 }
939
940 for (int icomp=0; icomp<3; icomp++) {
941 FillPatch(lev, t_old[lev], *mf_tke, GetVecOfPtrs(vec_tke), tke_bc(), BdyVars::null, icomp, false, false);
942 FillPatch(lev, t_old[lev], *mf_gls, GetVecOfPtrs(vec_gls), tke_bc(), BdyVars::null, icomp, false, false);
943 }
944 for (int icomp=0; icomp<NAT; icomp++) {
945 FillPatch(lev, t_old[lev], *mf_Akt, GetVecOfPtrs(vec_Akt), zvel_bc(), BdyVars::null, icomp, false, false);
946 }
947 FillPatchNoBC(lev, t_old[lev], *mf_Akv, GetVecOfPtrs(vec_Akv), BdyVars::null);
948 FillPatchNoBC(lev, t_old[lev], *mf_Akp, GetVecOfPtrs(vec_Akp), BdyVars::null);
949 FillPatchNoBC(lev, t_old[lev], *mf_Akk, GetVecOfPtrs(vec_Akk), BdyVars::null);
950}
two
constexpr amrex::Real two
Definition REMORA_Constants.H:8
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
vonKar
constexpr amrex::Real vonKar
Definition REMORA_Constants.H:38
GLS_StabilityType::Canuto_B
@ Canuto_B
GLS_StabilityType::Canuto_A
@ Canuto_A
GLS_StabilityType::Galperin
@ Galperin
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:30
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:10
NAT
#define NAT
Definition REMORA_IndexDefines.H:17
setVal
mf_h setVal(geomdata.ProbHi(2))
if
if(traditional)
Definition REMORA_InitAnalyticBathymetry_DogboneAnalytic.H:11
cff1
Real cff1
Definition REMORA_InitAnalyticBiology_BioToy.H:2
cff2
Real cff2
Definition REMORA_InitAnalyticBiology_BioToy.H:3
dhi
const auto dhi
Definition REMORA_InitAnalyticMasks_DogboneAnalytic.H:29
dlo
const auto dlo
Definition REMORA_InitAnalyticMasks_DogboneAnalytic.H:28
geomdata
auto geomdata
Definition REMORA_InitAnalyticProb_Advection.H:12
REMORA::zvel_bc
int zvel_bc() const noexcept
Definition REMORA.H:1487
REMORA::xvel_bc
int xvel_bc() const noexcept
Definition REMORA.H:1485
REMORA::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1772
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:603
REMORA::vec_ZoBot
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ZoBot
Bottom roughness length [m], defined at rho points.
Definition REMORA.H:537
REMORA::tke_bc
int tke_bc() const noexcept
Definition REMORA.H:1491
REMORA::vec_tke
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:667
REMORA::vec_gls
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:669
REMORA::vec_sustr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_sustr
Surface stress in the u direction.
Definition REMORA.H:479
REMORA::yvel_bc
int yvel_bc() const noexcept
Definition REMORA.H:1486
REMORA::foextrap_bc
int foextrap_bc() const noexcept
Definition REMORA.H:1493
REMORA::vec_Lscale
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Lscale
Vertical mixing turbulent length scale.
Definition REMORA.H:671
REMORA::vec_Hz
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hz
Width of cells in the vertical (z-) direction (3D, Hz in ROMS)
Definition REMORA.H:424
REMORA::vec_Akt
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:444
REMORA::xvel_old
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:386
REMORA::gls_prestep
void gls_prestep(int lev, amrex::MultiFab *mf_gls, amrex::MultiFab *mf_tke, amrex::MultiFab &mf_W, amrex::MultiFab *mf_msku, amrex::MultiFab *mf_mskv, const int nstp, const int nnew, const int iic, const int ntfirst, const int N, const amrex::Real dt_lev)
Prestep for GLS calculation.
Definition REMORA_gls.cpp:20
REMORA::vec_bvf
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvf
Brunt-Vaisala frequency (3D)
Definition REMORA.H:654
REMORA::physbcs
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1754
REMORA::vec_svstr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_svstr
Surface stress in the v direction.
Definition REMORA.H:481
REMORA::vec_Huon
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Huon
u-volume flux (3D)
Definition REMORA.H:426
REMORA::yvel_old
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:388
REMORA::t_new
amrex::Vector< amrex::Real > t_new
new time at each level, in seconds since start_time
Definition REMORA.H:1744
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:2017
REMORA::vec_Akk
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:673
REMORA::FillPatchNoBC
void FillPatchNoBC(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=false, amrex::Vector< amrex::MultiFab * > const &mfs_crse_old={}, amrex::Vector< amrex::MultiFab * > const &mfs_crse_new={})
Fill a new MultiFab by copying in phi from valid region and filling ghost cells without applying boun...
Definition REMORA_FillPatch.cpp:211
REMORA::FillPatch
void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=false, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=zero, const amrex::MultiFab &mf_calc=amrex::MultiFab(), amrex::Vector< amrex::MultiFab * > const &mfs_crse_old={}, amrex::Vector< amrex::MultiFab * > const &mfs_crse_new={})
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
Definition REMORA_FillPatch.cpp:30
REMORA::cons_old
amrex::Vector< amrex::MultiFab * > cons_old
multilevel data container for last step's scalar data: temperature, salinity, passive tracer
Definition REMORA.H:384
REMORA::vec_bustr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bustr
Bottom stress in the u direction.
Definition REMORA.H:540
REMORA::vec_bvstr
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvstr
Bottom stress in the v direction.
Definition REMORA.H:542
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:605
REMORA::vec_Akv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:442
REMORA::t_old
amrex::Vector< amrex::Real > t_old
old time at each level, in seconds since start_time
Definition REMORA.H:1746
REMORA::vec_z_w
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_w
z coordinates at w points (faces between z-cells)
Definition REMORA.H:456
REMORA::vec_Hvom
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hvom
v-volume flux (3D)
Definition REMORA.H:428
REMORA::vec_Akp
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:675
REMORA::gls_corrector
void gls_corrector(int lev, amrex::MultiFab *mf_gls, amrex::MultiFab *mf_tke, amrex::MultiFab &mf_W, amrex::MultiFab *mf_Akv, amrex::MultiFab *mf_Akt, amrex::MultiFab *mf_Akk, amrex::MultiFab *mf_Akp, amrex::MultiFab *mf_mskr, amrex::MultiFab *mf_msku, amrex::MultiFab *mf_mskv, const int nstp, const int nnew, const int N, const amrex::Real dt_lev)
Corrector step for GLS calculation.
Definition REMORA_gls.cpp:247
BdyVars::null
static constexpr int null
Definition REMORA_IndexDefines.H:69
SolverChoice::Akt_bak
amrex::Vector< amrex::Real > Akt_bak
Definition REMORA_DataStruct.H:1081
SolverChoice::gls_stability_type
GLS_StabilityType gls_stability_type
Definition REMORA_DataStruct.H:979
SolverChoice::Akv_bak
amrex::Real Akv_bak
Definition REMORA_DataStruct.H:1076
SolverChoice::gls_sigp
amrex::Real gls_sigp
Definition REMORA_DataStruct.H:1071
SolverChoice::gls_sigk
amrex::Real gls_sigk
Definition REMORA_DataStruct.H:1070
SolverChoice::gls_cmu0
amrex::Real gls_cmu0
Definition REMORA_DataStruct.H:1065
SolverChoice::Akk_bak
amrex::Real Akk_bak
Definition REMORA_DataStruct.H:1074
SolverChoice::gls_c3m
amrex::Real gls_c3m
Definition REMORA_DataStruct.H:1068
SolverChoice::gls_Gh0
amrex::Real gls_Gh0
Definition REMORA_DataStruct.H:1084
SolverChoice::gls_Ghmin
amrex::Real gls_Ghmin
Definition REMORA_DataStruct.H:1086
SolverChoice::gls_Kmin
amrex::Real gls_Kmin
Definition REMORA_DataStruct.H:1058
SolverChoice::gls_c3p
amrex::Real gls_c3p
Definition REMORA_DataStruct.H:1069
SolverChoice::gls_Ghcri
amrex::Real gls_Ghcri
Definition REMORA_DataStruct.H:1085
SolverChoice::gls_Pmin
amrex::Real gls_Pmin
Definition REMORA_DataStruct.H:1059
SolverChoice::Akp_bak
amrex::Real Akp_bak
Definition REMORA_DataStruct.H:1075