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