REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_uv3dmix.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] xbx x-nodal tile box
7 * @param[in ] ybx y-nodal tile box
8 * @param[inout] u u-direction velocity
9 * @param[inout] v v-direction velocity
10 * @param[in ] uold u-direction velocity at last step
11 * @param[in ] vold v-direction velocity at last step
12 * @param[inout] rufrc forcing term for u-velocity RHS
13 * @param[inout] rvfrc forcing term for v-velocity RHS
14 * @param[in ] visc2_p harmonic viscosity on psi points
15 * @param[in ] visc2_r harmonic viscosity on rho points
16 * @param[in ] Hz vertical cell height
17 * @param[in ] pm 1/dx
18 * @param[in ] pn 1/dy
19 * @param[in ] mskp land-sea mask on psi-points
20 * @param[in ] nrhs index of RHS component
21 * @param[in ] nnew index of time step to update
22 * @param[in ] dt_lev time step at this level
23 */
24void
25REMORA::uv3dmix (const Box& xbx, const Box& ybx,
26 const Array4<Real >& u,
27 const Array4<Real >& v,
28 const Array4<Real const>& uold,
29 const Array4<Real const>& vold,
30 const Array4<Real >& rufrc,
31 const Array4<Real >& rvfrc,
32 const Array4<Real const>& visc2_p,
33 const Array4<Real const>& visc2_r,
34 const Array4<Real const>& Hz,
35 const Array4<Real const>& pm,
36 const Array4<Real const>& pn,
37 const Array4<Real const>& mskp,
38 int nrhs, int nnew,
39 const Real dt_lev)
40{
41 //-----------------------------------------------------------------------
42 // Add in harmonic viscosity s terms.
43 //-----------------------------------------------------------------------
44
45 FArrayBox fab_UFx(growLo(xbx,0,1),1,amrex::The_Async_Arena()); fab_UFx.template setVal<RunOn::Device>(0.);
46 FArrayBox fab_UFe(growHi(xbx,1,1),1,amrex::The_Async_Arena()); fab_UFe.template setVal<RunOn::Device>(0.);
47 FArrayBox fab_VFe(growLo(ybx,1,1),1,amrex::The_Async_Arena()); fab_VFe.template setVal<RunOn::Device>(0.);
48 FArrayBox fab_VFx(growHi(ybx,0,1),1,amrex::The_Async_Arena()); fab_VFx.template setVal<RunOn::Device>(0.);
49
50 auto UFx=fab_UFx.array();
51 auto UFe=fab_UFe.array();
52 auto VFx=fab_VFx.array();
53 auto VFe=fab_VFe.array();
54
55 auto N = xbx.hiVect()[2] - ybx.loVect()[2];
56
57 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
58 //
59 // xbx is the x-face-centered box on which we update u [ 0:nx , 0:ny-1] (1,0,0) x-faces
60 // to do so requires UFx on [-1:nx , 0:ny-1] (0,0,0) cc
61 // which requires uold on [-1:nx+1, 0:ny-1] (1,0,0) x-faces
62 // and requires vold on [-1:nx , 0:ny ] (0,1,0) y-faces
63 // to do so requires UFe on [ 0:nx , 0:ny ] (1,1,0) xy-nodes
64 // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces
65 // and requires vold on [ -1:nx , 0:ny ] (0,1,0) y-faces
66
67 ParallelFor(growLo(xbx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
68 {
69 const Real cff = 0.5_rt*Hz(i,j,k) * ( pm(i,j,0) / pn(i,j,0) *
70 ( (pn(i ,j,0) + pn(i+1,j,0)) * uold(i+1,j,k,nrhs)-
71 (pn(i-1,j,0) + pn(i ,j,0)) * uold(i ,j,k,nrhs) )-
72 pn(i,j,0) / pm(i,j,0) *
73 ( (pm(i,j ,0) + pm(i,j+1,0)) * vold(i,j+1,k,nrhs)-
74 (pm(i,j-1,0) + pm(i,j ,0)) * vold(i,j ,k,nrhs) ) );
75
76 Real on_r = 1.0_rt / pn(i,j,0);
77 UFx(i,j,k) = on_r * on_r * visc2_r(i,j,0) * cff;
78 });
79
80 ParallelFor(growHi(xbx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
81 {
82 const Real pmon_p = (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0)) /
83 (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0));
84 const Real pnom_p = (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)) /
85 (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0));
86 const Real cff = mskp(i,j,0) * 0.125_rt *
87 (Hz(i-1,j ,k) + Hz(i,j ,k)+ Hz(i-1,j-1,k) + Hz(i,j-1,k))*
88 (pmon_p*
89 ((pn(i ,j-1,0)+pn(i ,j,0))*vold(i ,j,k,nrhs)-
90 (pn(i-1,j-1,0)+pn(i-1,j,0))*vold(i-1,j,k,nrhs))+
91 pnom_p*
92 ((pm(i-1,j ,0)+pm(i,j ,0))*uold(i,j ,k,nrhs)-
93 (pm(i-1,j-1,0)+pm(i,j-1,0))*uold(i,j-1,k,nrhs)));
94
95 const Real om_p = 4.0_rt / (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0));
96 UFe(i,j,k) = om_p*om_p*visc2_p(i,j,0)*cff;
97 });
98
99 ParallelFor(makeSlab(xbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
100 {
101 for (int k=0; k<=N; k++) {
102 const Real cff=dt_lev*0.25_rt*(pm(i-1,j,0)+pm(i,j,0))*(pn(i-1,j,0)+pn(i,j,0));
103 const Real cff1=0.5_rt*(pn(i-1,j,0)+pn(i,j,0))*(UFx(i,j ,k)-UFx(i-1,j,k));
104 const Real cff2=0.5_rt*(pm(i-1,j,0)+pm(i,j,0))*(UFe(i,j+1,k)-UFe(i ,j,k));
105 const Real cff3=cff*(cff1+cff2);
106 u(i,j,k,nnew)=u(i,j,k,nnew)+cff3;
107 rufrc(i,j,0) += cff1+cff2;
108 }
109 });
110
111
112 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
113 //
114 // ybx is the y-face-centered box on which we update v [ 0:nx-1, 0:ny ] (1,0,0) x-faces
115 // to do so requires VFe on [ 0:nx-1,-1:ny ] (0,0,0) cc
116 // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces
117 // and requires vold on [ 0:nx-1,-1:ny+1] (0,1,0) y-faces
118 // to do so requires VFx on [ 0:nx , 0:ny ] (1,1,0) xy-nodes
119 // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces
120 // and requires vold on [-1:nx , 0:ny ] (0,1,0) y-faces
121
122 ParallelFor(growLo(ybx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
123 {
124 // cff depends on k, but UFx and VFe will only be affected by the last cell?
125 const Real cff = 0.5_rt*Hz(i,j,k) * (pm(i,j,0) / pn(i,j,0) *
126 ((pn(i, j,0) + pn(i+1,j,0)) * uold(i+1,j,k,nrhs)-
127 (pn(i-1,j,0) + pn(i, j,0)) * uold(i ,j,k,nrhs))-
128 pn(i,j,0) / pm(i,j,0) *
129 ((pm(i,j ,0)+pm(i,j+1,0))*vold(i,j+1,k,nrhs)-
130 (pm(i,j-1,0)+pm(i,j ,0))*vold(i,j ,k,nrhs)));
131
132 Real om_r = 1.0_rt / pm(i,j,0);
133 VFe(i,j,k) = om_r * om_r * visc2_r(i,j,0) * cff;
134 });
135
136 ParallelFor(growHi(ybx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
137 {
138 const Real pmon_p = (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0)) /
139 (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0));
140 const Real pnom_p = (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)) /
141 (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0));
142 const Real cff = mskp(i,j,0) * 0.125_rt * (Hz(i-1,j ,k)+Hz(i,j ,k)+
143 Hz(i-1,j-1,k)+Hz(i,j-1,k))*
144 (pmon_p*
145 ((pn(i ,j-1,0)+pn(i ,j,0))*vold(i ,j,k,nrhs)-
146 (pn(i-1,j-1,0)+pn(i-1,j,0))*vold(i-1,j,k,nrhs))+
147 pnom_p*
148 ((pm(i-1,j ,0)+pm(i,j ,0))*uold(i,j ,k,nrhs)-
149 (pm(i-1,j-1,0)+pm(i,j-1,0))*uold(i,j-1,k,nrhs)));
150
151 const Real on_p = 4.0_rt / (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0));
152 VFx(i,j,k) = on_p*on_p*visc2_p(i,j,0)*cff;
153 });
154
155 ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
156 {
157 for (int k=0; k<=N; k++) {
158 const Real cff=dt_lev*0.25_rt*(pm(i,j,0)+pm(i,j-1,0))*(pn(i,j,0)+pn(i,j-1,0));
159 const Real cff1=0.5_rt*(pn(i,j-1,0)+pn(i,j,0))*(VFx(i+1,j,k)-VFx(i,j ,k));
160 const Real cff2=0.5_rt*(pm(i,j-1,0)+pm(i,j,0))*(VFe(i ,j,k)-VFe(i,j-1,k));
161 const Real cff3=cff*(cff1-cff2);
162
163 v(i,j,k,nnew)=v(i,j,k,nnew)+cff3;
164 rvfrc(i,j,0) += cff1-cff2;
165 }
166 });
167
168}
REMORA::uv3dmix
void uv3dmix(const amrex::Box &xbx, const amrex::Box &ybx, const amrex::Array4< amrex::Real > &u, const amrex::Array4< amrex::Real > &v, const amrex::Array4< amrex::Real const > &uold, const amrex::Array4< amrex::Real const > &vold, const amrex::Array4< amrex::Real > &rufrc, const amrex::Array4< amrex::Real > &rvfrc, const amrex::Array4< amrex::Real const > &visc2_p, const amrex::Array4< amrex::Real const > &visc2_r, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &mskp, int nrhs, int nnew, const amrex::Real dt_lev)
Harmonic viscosity.
Definition REMORA_uv3dmix.cpp:25