REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_rhs_uv_3d.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 ] xbx Box for operations on x-velocity
8 * @param[in ] ybx Box for operations on y-velocity
9 * @param[in ] uold u-velocity at last time step
10 * @param[in ] vold v-velocity at last time step
11 * @param[inout] ru u-velocity RHS
12 * @param[inout] rv v-velocity RHS
13 * @param[inout] rufrc forcing for u-velocity RHS
14 * @param[inout] rvfrc forcing for v-velocity RHS
15 * @param[in ] sustr u-direction surface momentum flux
16 * @param[in ] svstr v-direction surface momentum flux
17 * @param[in ] bustr u-direction bottom stress
18 * @param[in ] bvstr v-direction bottom stress
19 * @param[in ] Huon u-volume flux
20 * @param[in ] Hvom v-volume flux
21 * @param[in ] pm 1/dx
22 * @param[in ] pn 1/dy
23 * @param[in ] W vertical velocity
24 * @param FC temporary
25 * @param[in ] nrhs index of component for RHS
26 * @param[in ] N number of vertical levels
27 */
28void
29REMORA::rhs_uv_3d (int lev,
30 const Box& xbx, const Box& ybx,
31 const Array4<Real const>& uold ,
32 const Array4<Real const>& vold,
33 const Array4<Real >& ru,
34 const Array4<Real >& rv,
35 const Array4<Real >& rufrc,
36 const Array4<Real >& rvfrc,
37 const Array4<Real const>& sustr,
38 const Array4<Real const>& svstr,
39 const Array4<Real const>& bustr,
40 const Array4<Real const>& bvstr,
41 const Array4<Real const>& Huon,
42 const Array4<Real const>& Hvom,
43 const Array4<Real const>& pm,
44 const Array4<Real const>& pn,
45 const Array4<Real const>& W ,
46 const Array4<Real >& FC,
47 int nrhs, int N)
48{
49 BL_PROFILE("REMORA::rhs_uv_3d()");
50 const Box& domain = geom[lev].Domain();
51 const auto dlo = amrex::lbound(domain);
52 const auto dhi = amrex::ubound(domain);
53
54 GeometryData const& geomdata = geom[0].data();
55 bool is_periodic_in_x = geomdata.isPeriodic(0);
56 bool is_periodic_in_y = geomdata.isPeriodic(1);
57
58 //
59 // Scratch space
60 //
61 FArrayBox fab_UFx(growLo(xbx,0,1),1,amrex::The_Async_Arena()); fab_UFx.template setVal<RunOn::Device>(0.);
62 FArrayBox fab_UFe(growHi(xbx,1,1),1,amrex::The_Async_Arena()); fab_UFe.template setVal<RunOn::Device>(0.);
63 FArrayBox fab_VFe(growLo(ybx,1,1),1,amrex::The_Async_Arena()); fab_VFe.template setVal<RunOn::Device>(0.);
64 FArrayBox fab_VFx(growHi(ybx,0,1),1,amrex::The_Async_Arena()); fab_VFx.template setVal<RunOn::Device>(0.);
65
66 auto UFx=fab_UFx.array();
67 auto UFe=fab_UFe.array();
68 auto VFx=fab_VFx.array();
69 auto VFe=fab_VFe.array();
70
71 auto uv_hadv_scheme = solverChoice.uv_Hadv_scheme;
72
73 //check this////////////
74 const Real Gadv = -0.25_rt;
75
76 if (uv_hadv_scheme == AdvectionScheme::upstream3) {
77 // *************************************************************
78 // UPDATING U
79 // *************************************************************
80
81 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
82 //
83 // xbx is the x-face-centered box on which we update u (with ru) [ 0:nx , 0:ny-1] (1,0,0) x-faces
84 // to do so requires UFx on [-1:nx , 0:ny-1] (0,0,0) cc
85 // which requires uold on [-2:nx+2, 0:ny-1] (1,0,0) x-faces
86 // and requires Huon on [-2:nx+2, 0:ny-1] (0,0,0) x-faces
87 // to do so requires UFe on [ 0:nx , 0:ny ] (1,1,0) xy-nodes
88 // which requires uold on [ 0:nx ,-2:ny+1] (1,0,0) x-faces
89 // and requires Hvom on [-2:nx+1, 0:ny-1] (0,1,0) y-faces
90
91 //
92 // Define UFx, the x-fluxes at cell centers for updating u
93 // (Note that grow arguments are (bx, dir, ng)
94 //
95 ParallelFor(growLo(xbx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
96 {
97 Real cff1 = uold(i,j,k,nrhs)+uold(i+1,j,k,nrhs);
98
99 Real uxx_i = uold(i-1,j,k,nrhs)-2.0_rt*uold(i ,j,k,nrhs)+uold(i+1,j,k,nrhs);
100 Real uxx_ip1 = uold(i ,j,k,nrhs)-2.0_rt*uold(i+1,j,k,nrhs)+uold(i+2,j,k,nrhs);
101 // Upwinding
102
103 Real Huxx_i = Huon(i-1,j,k)-2.0_rt*Huon(i ,j,k)+Huon(i+1,j,k);
104 Real Huxx_ip1 = Huon(i ,j,k)-2.0_rt*Huon(i+1,j,k)+Huon(i+2,j,k);
105
106 if (i == dlo.x && !is_periodic_in_x) {
107 uxx_i = uxx_ip1;
108 Huxx_i = Huxx_ip1;
109 }
110 else if (i == dhi.x && !is_periodic_in_x) {
111 uxx_ip1 = uxx_i;
112 Huxx_ip1 = Huxx_i;
113 }
114
115 Real cff = (cff1 > 0.0_rt) ? uxx_i : uxx_ip1;
116
117 Real Huon_avg = (Huon(i,j,k) + Huon(i+1,j,k));
118
119 UFx(i,j,k) = 0.25_rt*(cff1+Gadv*cff) * ( Huon_avg + 0.5_rt*Gadv*(Huxx_i + Huxx_ip1) );
120 });
121
122 //
123 // Define UFe, the y-fluxes at nodes for updating u
124 //
125 ParallelFor(growHi(xbx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
126 {
127 Real cff1 = uold(i,j,k,nrhs) + uold(i ,j-1,k,nrhs);
128 Real cff2 = Hvom(i,j,k) + Hvom(i-1,j ,k);
129
130 Real uee_jm1 = uold(i,j-2,k,nrhs) - 2.0_rt*uold(i,j-1,k,nrhs) + uold(i ,j,k,nrhs);
131 Real uee_j = uold(i,j-1,k,nrhs) - 2.0_rt*uold(i,j ,k,nrhs) + uold(i,j+1,k,nrhs);
132
133 if (j == dlo.y and !is_periodic_in_y) {
134 uee_jm1 = uee_j;
135 } else if (j == dhi.y+1 and !is_periodic_in_y) {
136 uee_j = uee_jm1;
137 }
138
139 // Upwinding
140 Real cff = (cff2 > 0.0_rt) ? uee_jm1 : uee_j;
141
142 Real Hvxx_i = Hvom(i-1,j,k)-2.0_rt*Hvom(i ,j,k)+Hvom(i+1,j,k);
143 Real Hvxx_im1 = Hvom(i-2,j,k)-2.0_rt*Hvom(i-1,j,k)+Hvom(i ,j,k);
144
145 UFe(i,j,k) = 0.25_rt * (cff1+Gadv*cff)* (cff2+Gadv*0.5_rt*(Hvxx_i + Hvxx_im1));
146 });
147 } else if (uv_hadv_scheme == AdvectionScheme::centered2) {
148 ParallelFor(growLo(xbx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
149 {
150 UFx(i,j,k) = 0.25_rt * (uold(i,j,k,nrhs) + uold(i+1,j,k,nrhs)) * (Huon(i,j,k)+Huon(i+1,j,k));
151 });
152 ParallelFor(growHi(xbx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
153 {
154 UFe(i,j,k) = 0.25_rt * (uold(i,j-1,k,nrhs) + uold(i,j,k,nrhs)) * (Hvom(i-1,j,k)+Hvom(i,j,k));
155 });
156 }
157
158 //
159 // Define the RHS for u by differencing fluxes
160 //
161 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
162 {
163 ru(i,j,k,nrhs) -= ( (UFx(i,j,k)-UFx(i-1,j,k)) + (UFe(i,j+1,k)-UFe(i ,j,k)) );
164 });
165
166 if (uv_hadv_scheme == AdvectionScheme::upstream3) {
167 ParallelFor(surroundingNodes(xbx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
168 {
169 //-----------------------------------------------------------------------
170 // Add in vertical advection.
171 //-----------------------------------------------------------------------
172 Real cff1=9.0_rt/16.0_rt;
173 Real cff2=1.0_rt/16.0_rt;
174
175 if (k>1 && k<=N-1)
176 {
177 FC(i,j,k)=( cff1*(uold(i ,j,k-1,nrhs)+ uold(i,j,k ,nrhs))
178 -cff2*(uold(i ,j,k-2,nrhs)+ uold(i,j,k+1,nrhs)) )*
179 ( cff1*( W(i ,j,k)+ W(i-1,j,k))
180 -cff2*( W(i+1,j,k)+ W(i-2,j,k)) );
181 }
182 else // this needs to be split up so that the following can be concurrent
183 {
184 FC(i,j,N+1)=0.0_rt;
185
186 FC(i,j,N)=( cff1*(uold(i ,j,N-1,nrhs)+ uold(i,j,N ,nrhs))
187 -cff2*(uold(i ,j,N-2,nrhs)+ uold(i,j,N ,nrhs)) )*
188 ( cff1*( W(i ,j,N)+ W(i-1,j,N))
189 -cff2*( W(i+1,j,N)+ W(i-2,j,N)) );
190
191 FC(i,j,1)=( cff1*(uold(i ,j,0,nrhs)+ uold(i,j,1,nrhs))
192 -cff2*(uold(i ,j,0,nrhs)+ uold(i,j,2,nrhs)) )*
193 ( cff1*( W(i ,j,1)+ W(i-1,j,1))
194 -cff2*( W(i+1,j,1)+ W(i-2,j,1)) );
195 FC(i,j,0)=0.0_rt;
196 }
197 });
198 } else if (uv_hadv_scheme == AdvectionScheme::centered2) {
199 ParallelFor(surroundingNodes(xbx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
200 {
201 if (k>0 && k<=N) {
202 FC(i,j,k) = 0.25_rt * (uold(i,j,k-1,nrhs)+uold(i,j,k,nrhs)) * (W(i,j,k) + W(i-1,j,k));
203 } else {
204 FC(i,j,N+1)=0.0_rt;
205 FC(i,j,0)=0.0_rt;
206 }
207 });
208 }
209
210 ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
211 {
212 Real cff = FC(i,j,k+1)-FC(i,j,k);
213
214 ru(i,j,k,nrhs) -= cff;
215 });
216
217 Gpu::synchronize();
218
219 AMREX_ASSERT(xbx.smallEnd(2) == 0 && xbx.bigEnd(2) == N);
220 ParallelFor(makeSlab(xbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
221 {
222 for (int k = 0; k <= N; ++k)
223 {
224 rufrc(i,j,0) += ru(i,j,k,nrhs);
225
226 Real om_u = 2.0_rt / (pm(i-1,j,0)+pm(i,j,0));
227 Real on_u = 2.0_rt / (pn(i-1,j,0)+pn(i,j,0));
228 Real cff = om_u * on_u;
229
230 Real cff1 = (k == N) ? sustr(i,j,0)*cff : 0.0_rt;
231 Real cff2 = (k == 0) ? -bustr(i,j,0)*cff : 0.0_rt;
232
233 rufrc(i,j,0) += cff1+cff2;
234 }
235 });
236
237 if (uv_hadv_scheme == AdvectionScheme::upstream3) {
238 // *************************************************************
239 // UPDATING V
240 // *************************************************************
241
242 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
243 //
244 // ybx is the y-face-centered box on which we update v (with rv) [ 0:nx-1, 0:ny ] (1,0,0) y-faces
245 // to do so requires VFe on [ 0:nx-1,-1:ny ] (1,1,0) xy-nodes
246 // which requires vold on [ 0:nx-1,-2:ny+2] (1,0,0) y-faces
247 // and requires Hvom on [ 0:nx-1,-2:ny+2] (0,1,0) x-faces
248 // to do so requires VFx on [ 0:nx , 0:ny ] (0,0,0) cc
249 // which requires vold on [-2:nx+1, -:ny ] (1,0,0) y-faces
250 // and requires Hvom on [ 0:nx-1,-2:ny+1] (0,0,0) y-faces
251
252 // Grow ybx by one in high x-direction
253 ParallelFor(growHi(ybx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
254 {
255 Real cff1 = vold(i,j,k,nrhs) + vold(i-1,j ,k,nrhs);
256 Real cff2 = Huon(i,j,k) + Huon(i ,j-1,k);
257
258 Real vxx_im1 = vold(i-2,j,k,nrhs)-2.0_rt*vold(i-1,j,k,nrhs)+vold(i ,j,k,nrhs);
259 Real vxx_i = vold(i-1,j,k,nrhs)-2.0_rt*vold(i ,j,k,nrhs)+vold(i+1,j,k,nrhs);
260
261 if (i == dlo.x and !is_periodic_in_x) {
262 vxx_im1 = vxx_i;
263 } else if (i == dhi.x+1 and !is_periodic_in_x) {
264 vxx_i = vxx_im1;
265 }
266
267 // Upwinding
268 Real cff = (cff2 > 0.0_rt) ? vxx_im1 : vxx_i;
269
270
271 Real Huee_j = Huon(i,j-1,k)-2.0_rt*Huon(i,j ,k)+Huon(i,j+1,k);
272 Real Huee_jm1 = Huon(i,j-2,k)-2.0_rt*Huon(i,j-1,k)+Huon(i,j ,k);
273
274 VFx(i,j,k) = 0.25_rt*(cff1+Gadv*cff)* (cff2+Gadv*0.5_rt*(Huee_j + Huee_jm1));
275 });
276
277 // Grow ybx by one in low y-direction
278 ParallelFor(growLo(ybx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
279 {
280 Real cff1=vold(i,j,k,nrhs)+vold(i,j+1,k,nrhs);
281
282 // Upwinding
283 Real vee_j = vold(i,j-1,k,nrhs)-2.0_rt*vold(i,j ,k,nrhs)+ vold(i,j+1,k,nrhs);
284 Real vee_jp1 = vold(i,j ,k,nrhs)-2.0_rt*vold(i,j+1,k,nrhs)+ vold(i,j+2,k,nrhs);
285
286 Real Hvee_j = Hvom(i,j-1,k)-2.0_rt*Hvom(i,j ,k)+Hvom(i,j+1,k);
287 Real Hvee_jp1 = Hvom(i,j ,k)-2.0_rt*Hvom(i,j+1,k)+Hvom(i,j+2,k);
288
289 if (j == dlo.y and !is_periodic_in_y) {
290 vee_j = vee_jp1;
291 Hvee_j = Hvee_jp1;
292 }
293 else if (j == dhi.y and !is_periodic_in_y) {
294 vee_jp1 = vee_j;
295 Hvee_jp1 = Hvee_j;
296 }
297 Real cff = (cff1 > 0.0_rt) ? vee_j : vee_jp1;
298
299 VFe(i,j,k) = 0.25_rt * (cff1+Gadv*cff) * ( Hvom(i,j ,k)+ Hvom(i,j+1,k) + 0.5_rt * Gadv * (Hvee_j + Hvee_jp1) );
300 });
301 } else if (uv_hadv_scheme == AdvectionScheme::centered2) {
302 ParallelFor(growHi(ybx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
303 {
304 VFx(i,j,k) = 0.25_rt * (vold(i-1,j,k,nrhs) + vold(i,j,k,nrhs)) * (Huon(i,j-1,k)+Huon(i,j,k));
305 });
306 ParallelFor(growLo(ybx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
307 {
308 VFe(i,j,k) = 0.25_rt * (vold(i,j,k,nrhs) + vold(i,j+1,k,nrhs)) * (Hvom(i,j,k)+Hvom(i,j+1,k));
309 });
310 }
311
312 ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
313 {
314 rv(i,j,k,nrhs) -= ( (VFx(i+1,j,k)-VFx(i,j,k)) + (VFe(i,j,k)-VFe(i,j-1,k)) );
315 });
316
317 if (uv_hadv_scheme == AdvectionScheme::upstream3) {
318 ParallelFor(surroundingNodes(ybx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
319 {
320 Real cff1=9.0_rt/16.0_rt;
321 Real cff2=1.0_rt/16.0_rt;
322
323 if (k>1 && k<=N-1)
324 {
325 FC(i,j,k)=( cff1*(vold(i,j,k-1,nrhs)+ vold(i,j,k ,nrhs))
326 -cff2*(vold(i,j,k-2,nrhs)+ vold(i,j,k+1,nrhs)) )*
327 ( cff1*(W(i,j ,k)+ W(i,j-1,k))
328 -cff2*(W(i,j+1,k)+ W(i,j-2,k)) );
329 }
330 else // this needs to be split up so that the following can be concurrent
331 {
332 FC(i,j,N+1)=0.0_rt;
333 FC(i,j,N)=( cff1*(vold(i,j,N-1,nrhs)+ vold(i,j,N ,nrhs))
334 -cff2*(vold(i,j,N-2,nrhs)+ vold(i,j,N ,nrhs)) )*
335 ( cff1*(W(i,j ,N)+ W(i,j-1,N))
336 -cff2*(W(i,j+1,N)+ W(i,j-2,N)) );
337 FC(i,j,1)=( cff1*(vold(i,j,0,nrhs)+ vold(i,j,1,nrhs))
338 -cff2*(vold(i,j,0,nrhs)+ vold(i,j,2,nrhs)) )*
339 ( cff1*(W(i,j ,1)+ W(i,j-1,1))
340 -cff2*(W(i,j+1,1)+ W(i,j-2,1)) );
341 FC(i,j,0)=0.0_rt;
342 // FC(i,0,-1)=0.0_rt;
343 }
344 });
345 } else if (uv_hadv_scheme == AdvectionScheme::centered2) {
346 ParallelFor(surroundingNodes(ybx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
347 {
348 if (k>0 && k<=N) {
349 FC(i,j,k) = 0.25_rt * (vold(i,j,k-1,nrhs)+vold(i,j,k,nrhs)) * (W(i,j,k) + W(i,j-1,k));
350 } else {
351 FC(i,j,N+1)=0.0_rt;
352 FC(i,j,0)=0.0_rt;
353 }
354 });
355 }
356 Gpu::synchronize();
357
358 ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
359 {
360 Real cff = FC(i,j,k+1)-FC(i,j,k);
361
362 rv(i,j,k,nrhs) -= cff;
363 });
364
365 Gpu::synchronize();
366
367 ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
368 {
369 for (int k = 0; k <= N; ++k)
370 {
371 rvfrc(i,j,0) += rv(i,j,k,nrhs);
372
373 Real om_v = 2.0_rt / (pm(i,j-1,0)+pm(i,j,0));
374 Real on_v = 2.0_rt / (pn(i,j-1,0)+pn(i,j,0));
375 Real cff = om_v * on_v;
376
377 Real cff1 = (k == N) ? svstr(i,j,0)*cff : 0.0_rt;
378 Real cff2 = (k == 0) ? -bvstr(i,j,0)*cff : 0.0_rt;
379
380 rvfrc(i,j,0) += cff1+cff2;
381 }
382 });
383}
AdvectionScheme::centered2
@ centered2
AdvectionScheme::upstream3
@ upstream3
REMORA::rhs_uv_3d
void rhs_uv_3d(int lev, const amrex::Box &xbx, const amrex::Box &ybx, const amrex::Array4< amrex::Real const > &uold, const amrex::Array4< amrex::Real const > &vold, const amrex::Array4< amrex::Real > &ru, const amrex::Array4< amrex::Real > &rv, const amrex::Array4< amrex::Real > &rufrc, const amrex::Array4< amrex::Real > &rvfrc, const amrex::Array4< amrex::Real const > &sustr, const amrex::Array4< amrex::Real const > &svstr, const amrex::Array4< amrex::Real const > &bustr, const amrex::Array4< amrex::Real const > &bvstr, const amrex::Array4< amrex::Real const > &Huon, const amrex::Array4< amrex::Real const > &Hvom, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &W, const amrex::Array4< amrex::Real > &FC, int nrhs, int N)
RHS terms for 3D momentum.
Definition REMORA_rhs_uv_3d.cpp:29
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
SolverChoice::uv_Hadv_scheme
AdvectionScheme uv_Hadv_scheme
Definition REMORA_DataStruct.H:572