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 BL_PROFILE("REMORA::uv3dmix()");
42 //-----------------------------------------------------------------------
43 // Add in harmonic viscosity s terms.
44 //-----------------------------------------------------------------------
45
46 int ncomp = 0;
47 int UFx_comp = ncomp++;
48 int UFe_comp = ncomp++;
49 int VFe_comp = ncomp++;
50 int VFx_comp = ncomp++;
51
52 Box bx = enclosedCells(xbx);
53 FArrayBox fab(grow(bx,IntVect(2,2,0)),ncomp,amrex::The_Async_Arena()); fab.template setVal<RunOn::Device>(0.);
54
55 auto UFx=fab.array(UFx_comp);
56 auto UFe=fab.array(UFe_comp);
57 auto VFx=fab.array(VFx_comp);
58 auto VFe=fab.array(VFe_comp);
59
60 auto N = xbx.hiVect()[2] - ybx.loVect()[2];
61
62 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
63 //
64 // xbx is the x-face-centered box on which we update u [ 0:nx , 0:ny-1] (1,0,0) x-faces
65 // to do so requires UFx on [-1:nx , 0:ny-1] (0,0,0) cc
66 // which requires uold on [-1:nx+1, 0:ny-1] (1,0,0) x-faces
67 // and requires vold on [-1:nx , 0:ny ] (0,1,0) y-faces
68 // to do so requires UFe on [ 0:nx , 0:ny ] (1,1,0) xy-nodes
69 // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces
70 // and requires vold on [ -1:nx , 0:ny ] (0,1,0) y-faces
71
72 ParallelFor(growLo(xbx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
73 {
74 const Real cff = 0.5_rt*Hz(i,j,k) * ( pm(i,j,0) / pn(i,j,0) *
75 ( (pn(i ,j,0) + pn(i+1,j,0)) * uold(i+1,j,k,nrhs)-
76 (pn(i-1,j,0) + pn(i ,j,0)) * uold(i ,j,k,nrhs) )-
77 pn(i,j,0) / pm(i,j,0) *
78 ( (pm(i,j ,0) + pm(i,j+1,0)) * vold(i,j+1,k,nrhs)-
79 (pm(i,j-1,0) + pm(i,j ,0)) * vold(i,j ,k,nrhs) ) );
80
81 Real on_r = 1.0_rt / pn(i,j,0);
82 UFx(i,j,k) = on_r * on_r * visc2_r(i,j,0) * cff;
83 });
84
85 ParallelFor(growHi(xbx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
86 {
87 const Real pmon_p = (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0)) /
88 (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0));
89 const Real pnom_p = (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)) /
90 (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0));
91 const Real cff = mskp(i,j,0) * 0.125_rt *
92 (Hz(i-1,j ,k) + Hz(i,j ,k)+ Hz(i-1,j-1,k) + Hz(i,j-1,k))*
93 (pmon_p*
94 ((pn(i ,j-1,0)+pn(i ,j,0))*vold(i ,j,k,nrhs)-
95 (pn(i-1,j-1,0)+pn(i-1,j,0))*vold(i-1,j,k,nrhs))+
96 pnom_p*
97 ((pm(i-1,j ,0)+pm(i,j ,0))*uold(i,j ,k,nrhs)-
98 (pm(i-1,j-1,0)+pm(i,j-1,0))*uold(i,j-1,k,nrhs)));
99
100 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));
101 UFe(i,j,k) = om_p*om_p*visc2_p(i,j,0)*cff;
102 });
103
104 ParallelFor(makeSlab(xbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
105 {
106 for (int k=0; k<=N; k++) {
107 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));
108 const Real cff1=0.5_rt*(pn(i-1,j,0)+pn(i,j,0))*(UFx(i,j ,k)-UFx(i-1,j,k));
109 const Real cff2=0.5_rt*(pm(i-1,j,0)+pm(i,j,0))*(UFe(i,j+1,k)-UFe(i ,j,k));
110 const Real cff3=cff*(cff1+cff2);
111 u(i,j,k,nnew)=u(i,j,k,nnew)+cff3;
112 rufrc(i,j,0) += cff1+cff2;
113 }
114 });
115
116
117 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
118 //
119 // ybx is the y-face-centered box on which we update v [ 0:nx-1, 0:ny ] (1,0,0) x-faces
120 // to do so requires VFe on [ 0:nx-1,-1:ny ] (0,0,0) cc
121 // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces
122 // and requires vold on [ 0:nx-1,-1:ny+1] (0,1,0) y-faces
123 // to do so requires VFx on [ 0:nx , 0:ny ] (1,1,0) xy-nodes
124 // which requires uold on [ 0:nx ,-1:ny ] (1,0,0) x-faces
125 // and requires vold on [-1:nx , 0:ny ] (0,1,0) y-faces
126
127 ParallelFor(growLo(ybx,1,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
128 {
129 // cff depends on k, but UFx and VFe will only be affected by the last cell?
130 const Real cff = 0.5_rt*Hz(i,j,k) * (pm(i,j,0) / pn(i,j,0) *
131 ((pn(i, j,0) + pn(i+1,j,0)) * uold(i+1,j,k,nrhs)-
132 (pn(i-1,j,0) + pn(i, j,0)) * uold(i ,j,k,nrhs))-
133 pn(i,j,0) / pm(i,j,0) *
134 ((pm(i,j ,0)+pm(i,j+1,0))*vold(i,j+1,k,nrhs)-
135 (pm(i,j-1,0)+pm(i,j ,0))*vold(i,j ,k,nrhs)));
136
137 Real om_r = 1.0_rt / pm(i,j,0);
138 VFe(i,j,k) = om_r * om_r * visc2_r(i,j,0) * cff;
139 });
140
141 ParallelFor(growHi(ybx,0,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
142 {
143 const Real pmon_p = (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0)) /
144 (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0));
145 const Real pnom_p = (pn(i-1,j-1,0)+pn(i-1,j,0)+pn(i,j-1,0)+pn(i,j,0)) /
146 (pm(i-1,j-1,0)+pm(i-1,j,0)+pm(i,j-1,0)+pm(i,j,0));
147 const Real cff = mskp(i,j,0) * 0.125_rt * (Hz(i-1,j ,k)+Hz(i,j ,k)+
148 Hz(i-1,j-1,k)+Hz(i,j-1,k))*
149 (pmon_p*
150 ((pn(i ,j-1,0)+pn(i ,j,0))*vold(i ,j,k,nrhs)-
151 (pn(i-1,j-1,0)+pn(i-1,j,0))*vold(i-1,j,k,nrhs))+
152 pnom_p*
153 ((pm(i-1,j ,0)+pm(i,j ,0))*uold(i,j ,k,nrhs)-
154 (pm(i-1,j-1,0)+pm(i,j-1,0))*uold(i,j-1,k,nrhs)));
155
156 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));
157 VFx(i,j,k) = on_p*on_p*visc2_p(i,j,0)*cff;
158 });
159
160 ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int )
161 {
162 for (int k=0; k<=N; k++) {
163 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));
164 const Real cff1=0.5_rt*(pn(i,j-1,0)+pn(i,j,0))*(VFx(i+1,j,k)-VFx(i,j ,k));
165 const Real cff2=0.5_rt*(pm(i,j-1,0)+pm(i,j,0))*(VFe(i ,j,k)-VFe(i,j-1,k));
166 const Real cff3=cff*(cff1-cff2);
167
168 v(i,j,k,nnew)=v(i,j,k,nnew)+cff3;
169 rvfrc(i,j,0) += cff1-cff2;
170 }
171 });
172}
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