REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_prestep_t_advection.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[in ] tbx tile box
8 * @param[in ] gbx grown tile box
9 * @param[in ] tempold scalar at last time
10 * @param[in ] tempcache cached current time step's scalar value
11 * @param[in ] Hz vertical cell height
12 * @param[in ] Huon u-volume flux
13 * @param[in ] Hvom v-volume flux
14 * @param[in ] Akv vertical viscosity coefficient
15 * @param[inout] W vertical velocity
16 * @param DC temporary
17 * @param FC temporary
18 * @param[ out] tempstore scratch space for calculations on scalars
19 * @param[in ] z_w z coordinates at w points
20 * @param[in ] h bathymetry
21 * @param[in ] pm 1/dx
22 * @param[in ] pn 1/dy
23 * @param[in ] msku land-sea mask on u-points
24 * @param[in ] mskv land-sea mask on v-points
25 * @param[in ] river_pos river positions
26 * @param[in ] river_source river source data to add, if using
27 * @param[in ] iic which time step we're on
28 * @param[in ] ntfirst what is the first time step?
29 * @param[in ] nrhs index of RHS component
30 * @param[in ] N number of vertical levels
31 * @param[in ] dt_lev time step at this level
32 */
33
34void
35REMORA::prestep_t_advection (int lev, const Box& tbx, const Box& gbx,
36 const Array4<Real >& tempold,
37 const Array4<Real >& tempcache,
38 const Array4<Real >& Hz,
39 const Array4<Real >& Huon,
40 const Array4<Real >& Hvom,
41 const Array4<Real >& W ,
42 const Array4<Real >& DC,
43 const Array4<Real >& FC ,
44 const Array4<Real >& tempstore,
45 const Array4<Real const>& z_w,
46 const Array4<Real const>& h,
47 const Array4<Real const>& pm,
48 const Array4<Real const>& pn,
49 const Array4<Real const>& msku,
50 const Array4<Real const>& mskv,
51 const Array4<int const>& river_pos,
52 const Array4<Real const>& river_source,
53 int iic, int ntfirst, int nrhs, int N,
54 Real dt_lev)
55{
56 BL_PROFILE("REMORA::prestep_t_advection()");
57 const Box& domain = geom[lev].Domain();
58 const auto dlo = amrex::lbound(domain);
59 const auto dhi = amrex::ubound(domain);
60
61 GeometryData const& geomdata = geom[0].data();
62 bool is_periodic_in_x = geomdata.isPeriodic(0);
63 bool is_periodic_in_y = geomdata.isPeriodic(1);
64
65 //copy the tilebox
66 Box gbx1 = tbx;
67 Box gbx2 = tbx;
68 Box tbxp1 = tbx;
69 Box tbxp2 = tbx;
70
71 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
72 tbxp2.grow(IntVect(NGROW,NGROW,0));
73
74 int ncomp = 0;
75 int FX_comp = ncomp++;
76 int FE_comp = ncomp++;
77 int curv_comp = ncomp++;
78 int grad_comp = ncomp++;
79
80 FArrayBox fab(tbxp2,ncomp,amrex::The_Async_Arena());
81 fab.template setVal<RunOn::Device>(0.0_rt);
82
83 auto FX=fab.array(FX_comp);
84 auto FE=fab.array(FE_comp);
85 auto curv=fab.array(curv_comp);
86 auto grad=fab.array(grad_comp);
87
88 Box ubx = surroundingNodes(tbx,0);
89 Box vbx = surroundingNodes(tbx,1);
90
91 Box utbxp1 = surroundingNodes(tbxp1,0);
92 Box vtbxp1 = surroundingNodes(tbxp1,1);
93
94 Box gbx3uneven_init(IntVect(AMREX_D_DECL(tbx.smallEnd(0)-3,tbx.smallEnd(1)-3,tbx.smallEnd(2))),
95 IntVect(AMREX_D_DECL(tbx.bigEnd(0)+2,tbx.bigEnd(1)+2,tbx.bigEnd(2))));
96 BoxArray ba_gbx3uneven = intersect(BoxArray(gbx3uneven_init), gbx);
97 AMREX_ASSERT((ba_gbx3uneven.size() == 1));
98 Box gbx3uneven = ba_gbx3uneven[0];
99
100 gbx2.grow(IntVect(NGROW,NGROW,0));
101 BoxArray ba_gbx2 = intersect(BoxArray(gbx2), gbx);
102 AMREX_ASSERT((ba_gbx2.size() == 1));
103 gbx2 = ba_gbx2[0];
104
105 gbx1.grow(IntVect(NGROW-1,NGROW-1,0));
106 BoxArray ba_gbx1 = intersect(BoxArray(gbx1), gbx);
107 AMREX_ASSERT((ba_gbx1.size() == 1));
108 gbx1 = ba_gbx1[0];
109
110 //------------------------------------------------------------------------
111 // Vertically integrate horizontal mass flux divergence.
112 //------------------------------------------------------------------------
113 //
114 //Should really use gbx3uneven
115 Box gbx3unevenD = gbx3uneven;
116 gbx3unevenD.makeSlab(2,0);
117 Box gbx1D = gbx1;
118 gbx1D.makeSlab(2,0);
119
120 // We used to set W(i,j,0) = 0.0, but this should have already been set before passing into the function
121 ParallelFor(gbx1D, N+1,
122 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
123 {
124 // Starting with zero vertical velocity at the bottom, integrate
125 // from the bottom (k=0) to the free-surface (k=N). The w(:,:,N(ng))
126 // contains the vertical velocity at the free-surface, d(zeta)/d(t).
127 // Notice that barotropic mass flux divergence is not used directly.
128 //
129 int k = kk + 1;
130 W(i,j,k) = W(i,j,k-1) - (Huon(i+1,j,k-1)-Huon(i,j,k-1)) - (Hvom(i,j+1,k-1)-Hvom(i,j,k-1));
131 });
132 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
133 {
134 W(i,j,N+1)=W(i,j,N+1)/(z_w(i,j,N+1)+h(i,j,0,0)); // wrk_i
135 });
136 ParallelFor(gbx1D, N,
137 [=] AMREX_GPU_DEVICE (int i, int j, int , int kk)
138 {
139 int k = kk + 1;
140 W(i,j,k) = W(i,j,k)- W(i,j,N+1)*(z_w(i,j,k)+h(i,j,0,0));
141 });
142 ParallelFor(gbx1D, [=] AMREX_GPU_DEVICE (int i, int j, int )
143 {
144 W(i,j,N+1) = 0.0_rt;
145 });
146
147 //Use FC and DC as intermediate arrays for FX and FE
148 //First pass do centered 2d terms
149 if (solverChoice.flat_bathymetry) {
150 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
151 {
152 FX(i,j,k)=Box(tempold).contains(i-1,j,k) ? Huon(i,j,k)*
153 0.5_rt*(tempold(i-1,j,k)+tempold(i ,j,k)) : 1e34_rt;
154 });
155 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
156 {
157 FE(i,j,k)=Box(tempold).contains(i,j-1,k) ? Hvom(i,j,k)*
158 0.5_rt*(tempold(i,j-1,k)+tempold(i,j,k)) : 1e34_rt;
159 });
160
161 } else {
162
163 ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
164 {
165 //should be t index 3
166 FX(i,j,k)=(tempold(i,j,k,nrhs)-tempold(i-1,j,k,nrhs)) * msku(i,j,0);
167 });
168
169 Box utbxp1_slab_lo = makeSlab(utbxp1,0,dlo.x-1) & utbxp1;
170 Box utbxp1_slab_hi = makeSlab(utbxp1,0,dhi.x+1) & utbxp1;
171 if (utbxp1_slab_lo.ok() && !is_periodic_in_x) {
172 ParallelFor(utbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
173 {
174 FX(i,j,k) = FX(i+1,j,k);
175 });
176 }
177 if (utbxp1_slab_hi.ok() && !is_periodic_in_x) {
178 ParallelFor(utbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
179 {
180 FX(i+1,j,k) = FX(i,j,k);
181 });
182 }
183
184 ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
185 {
186 //should be t index 3
187 FE(i,j,k)=(tempold(i,j,k,nrhs)-tempold(i,j-1,k,nrhs)) * mskv(i,j,0);
188 });
189
190 Box vtbxp1_slab_lo = makeSlab(vtbxp1,1,dlo.y-1) & vtbxp1;
191 Box vtbxp1_slab_hi = makeSlab(vtbxp1,1,dhi.y+1) & vtbxp1;
192 if (vtbxp1_slab_lo.ok() && !is_periodic_in_y) {
193 ParallelFor(vtbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
194 {
195 FE(i,j,k) = FE(i,j+1,k);
196 });
197 }
198 if (vtbxp1_slab_hi.ok() && !is_periodic_in_y) {
199 ParallelFor(vtbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
200 {
201 FE(i,j+1,k) = FE(i,j,k);
202 });
203 }
204
205 Real cffa=1.0_rt/6.0_rt;
206 Real cffb=1.0_rt/3.0_rt;
207 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3)
208 {
209 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
210 {
211 //Upstream3
212 curv(i,j,k)=-FX(i,j,k)+FX(i+1,j,k);
213 });
214
215 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
216 {
217 Real max_Huon = std::max(Huon(i,j,k),0.0_rt);
218 Real min_Huon = std::min(Huon(i,j,k),0.0_rt);
219
220 FX(i,j,k)=Huon(i,j,k)*0.5_rt*(tempold(i,j,k)+tempold(i-1,j,k))-
221 cffa*(curv(i,j,k)*min_Huon+ curv(i-1,j,k)*max_Huon);
222 });
223
224 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
225
226 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
227 {
228 //Centered4
229 grad(i,j,k)=0.5_rt*(FX(i,j,k)+FX(i+1,j,k));
230 });
231
232 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
233 {
234 FX(i,j,k)=Huon(i,j,k)*0.5_rt*(tempold(i,j,k)+tempold(i-1,j,k)-
235 cffb*(grad(i,j,k)-grad(i-1,j,k)));
236 });
237 } else {
238 Error("Not a valid horizontal advection scheme");
239 }
240
241 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3)
242 {
243 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
244 {
245 curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k);
246 });
247
248 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
249 {
250 Real max_Hvom = std::max(Hvom(i,j,k),0.0_rt);
251 Real min_Hvom = std::min(Hvom(i,j,k),0.0_rt);
252
253 FE(i,j,k)=Hvom(i,j,k)*0.5_rt*(tempold(i,j,k)+tempold(i,j-1,k))-
254 cffa*(curv(i,j,k)*min_Hvom+ curv(i,j-1,k)*max_Hvom);
255 });
256
257 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
258
259 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
260 {
261 grad(i,j,k)=0.5_rt*(FE(i,j,k)+FE(i,j+1,k));
262 });
263
264 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
265 {
266 FE(i,j,k)=Hvom(i,j,k)*0.5_rt*(tempold(i,j,k)+tempold(i,j-1,k)-
267 cffb*(grad(i,j,k)- grad(i,j-1,k)));
268 });
269 } else {
270 Error("Not a valid horizontal advection scheme");
271 }
272 } // not flat
273
274 bool do_rivers_cons = (river_source.size() > 0);
275 if (solverChoice.do_rivers) {
276 int* river_direction_d = river_direction.data();
277 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
278 {
279 int iriver = river_pos(i,j,0);
280 if (iriver >= 0) {
281 if (river_direction_d[iriver] == 0) {
282 FX(i,j,k) = (!do_rivers_cons) ? 0.0_rt : Huon(i,j,k) * river_source(iriver,0,k);
283 } else {
284 FE(i,j,k) = (!do_rivers_cons) ? 0.0_rt : Hvom(i,j,k) * river_source(iriver,0,k);
285 }
286 }
287 });
288 }
289
290 //Intermediate tracer at 3
291 //
292 // Time-step horizontal advection (m Tunits).
293 //
294 Real cff1 = 0.0_rt, cff2 = 0.0_rt, cff;
295
296 Real GammaT = 1.0_rt/6.0_rt;
297
298 if (iic==ntfirst)
299 {
300 cff=0.5_rt*dt_lev;
301 cff1=1.0_rt;
302 cff2=0.0_rt;
303 } else {
304 cff=(1.0_rt-GammaT)*dt_lev;
305 cff1=0.5_rt+GammaT;
306 cff2=0.5_rt-GammaT;
307 }
308
309 ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
310 {
311 tempstore(i,j,k)=Hz(i,j,k)*( cff1 * tempold(i,j,k)+
312 cff2 * tempcache(i,j,k) )-
313 cff * pm(i,j,0)*pn(i,j,0) * (FX(i+1,j,k)-FX(i,j,k)+
314 FE(i,j+1,k)-FE(i,j,k));
315 });
316
317 Real c1=0.5_rt;
318 Real c2=7.0_rt/12.0_rt;
319 Real c3=1.0_rt/12.0_rt;
320
321 //
322 // Time-step vertical advection of tracers (Tunits). Impose artificial
323 // continuity equation.
324 //
325 ParallelFor(growHi(growLo(convert(tbx,IntVect(0,0,1)),2,-2),2,-1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
326 {
327 //-----------------------------------------------------------------------
328 // Add in vertical advection.
329 //-----------------------------------------------------------------------
330
331 FC(i,j,k)=( c2*(tempold(i ,j,k-1,nrhs)+ tempold(i,j,k ,nrhs))
332 -c3*(tempold(i ,j,k-2 ,nrhs)+ tempold(i,j,k+1,nrhs)) )*
333 ( W(i,j,k));
334 });
335 ParallelFor(makeSlab(tbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
336 {
337 FC(i,j,N+1)=0.0_rt;
338 FC(i,j,N) = ( c2*tempold(i,j,N-1,nrhs)+ c1*tempold(i,j,N,nrhs)-c3*tempold(i,j,N-2,nrhs) )
339 * W(i,j,N);
340 FC(i,j, 1) = ( c2*tempold(i,j, 1,nrhs)+ c1*tempold(i,j,0,nrhs)-c3*tempold(i,j,2,nrhs) )
341 * W(i,j,1);
342 FC(i,j, 0) = 0.0_rt;
343 });
344
345 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
346 {
347 DC(i,j,k)=1.0_rt/(Hz(i,j,k)-
348 cff*pm(i,j,0)*pn(i,j,0)*
349 (Huon(i+1,j,k)-Huon(i,j,k)+
350 Hvom(i,j+1,k)-Hvom(i,j,k)+
351 (W(i,j,k+1)-W(i,j,k))));
352 });
353
354 ParallelFor(tbx,
355 [=] AMREX_GPU_DEVICE (int i, int j, int k)
356 {
357 Real c_p = cff*pm(i,j,0)*pn(i,j,0);
358
359 Real c4 = FC(i,j,k+1)-FC(i,j,k);
360
361 tempstore(i,j,k) = DC(i,j,k)*(tempstore(i,j,k)-c_p*c4);
362 });
363}
AdvectionScheme::upstream3
@ upstream3
AdvectionScheme::centered4
@ centered4
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:13
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
REMORA::river_direction
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1096
REMORA::prestep_t_advection
void prestep_t_advection(int lev, const amrex::Box &tbx, const amrex::Box &gbx, const amrex::Array4< amrex::Real > &tempold, const amrex::Array4< amrex::Real > &tempcache, const amrex::Array4< amrex::Real > &Hz, const amrex::Array4< amrex::Real > &Huon, const amrex::Array4< amrex::Real > &Hvom, const amrex::Array4< amrex::Real > &W, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real > &sstore, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &h, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &msku, const amrex::Array4< amrex::Real const > &mskv, const amrex::Array4< int const > &river_pos, const amrex::Array4< amrex::Real const > &river_source, int iic, int ntfirst, int nrhs, int N, const amrex::Real dt_lev)
Prestep advection calculations for the tracers.
Definition REMORA_prestep_t_advection.cpp:35
SolverChoice::tracer_Hadv_scheme
AdvectionScheme tracer_Hadv_scheme
Definition REMORA_DataStruct.H:571