REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_rhs_uv_2d.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 ] ubar barotropic x-velocity
10 * @param[in ] vbar barotropic y-velocity
11 * @param[inout] rhs_ubar RHS for barotropic x-velocity
12 * @param[inout] rhs_vbar RHS for barotropic y-velocity
13 * @param[in ] DUon u-volume flux (barotropic)
14 * @param[in ] DVom v-volume flux (barotropic)
15 * @param[in ] krhs index of rhs component
16 */
17void
18REMORA::rhs_uv_2d (int lev, const Box& xbx, const Box& ybx,
19 const Array4<Real const>& ubar,
20 const Array4<Real const>& vbar,
21 const Array4<Real >& rhs_ubar ,
22 const Array4<Real >& rhs_vbar,
23 const Array4<Real const>& DUon,
24 const Array4<Real const>& DVom,
25 const int krhs)
26{
27 BL_PROFILE("REMORA::rhs_uv_2d()");
28 const Box& domain = geom[lev].Domain();
29 const auto dlo = amrex::lbound(domain);
30 const auto dhi = amrex::ubound(domain);
31
32 GeometryData const& geomdata = geom[0].data();
33 bool is_periodic_in_x = geomdata.isPeriodic(0);
34 bool is_periodic_in_y = geomdata.isPeriodic(1);
35
36 int ncompbc = 1;
37 Vector<BCRec> bcrs_x(ncompbc);
38 Vector<BCRec> bcrs_y(ncompbc);
39 amrex::setBC(xbx,domain,xvel_bc(),0,1,domain_bcs_type,bcrs_x);
40 amrex::setBC(ybx,domain,yvel_bc(),0,1,domain_bcs_type,bcrs_y);
41
42 //
43 // Scratch space
44 //
45
46 Box bx = enclosedCells(xbx);
47 int ncomp = 0;
48 int UFx_comp = ncomp++;
49 int UFe_comp = ncomp++;
50 int VFe_comp = ncomp++;
51 int VFx_comp = ncomp++;
52 FArrayBox fab(grow(bx,IntVect(2,2,0)),ncomp,amrex::The_Async_Arena());
53
54 auto UFx=fab.array(UFx_comp);
55 auto UFe=fab.array(UFe_comp);
56 auto VFx=fab.array(VFx_comp);
57 auto VFe=fab.array(VFe_comp);
58
59 auto uv_hadv_scheme = solverChoice.uv_Hadv_scheme;
60
61 if (uv_hadv_scheme == AdvectionScheme::upstream3) {
62 // *************************************************************
63 // UPDATING U
64 // *************************************************************
65
66 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
67 //
68 // xbx is the x-face-centered box on which we update ubar (with rhs_ubar) [ 0:nx , 0:ny-1] (1,0,0) x-faces
69 // to do so requires UFx on [-1:nx , 0:ny-1] (0,0,0) cc
70 // which requires ubar on [-2:nx+2, 0:ny-1] (1,0,0) x-faces
71 // and requires DUon on [-2:nx+2, 0:ny-1] (1,0,0) x-faces
72 // to do so requires UFe on [ 0:nx , 0:ny ] (1,1,0) xy-nodes
73 // which requires ubar on [ 0:nx ,-2:ny+1] (1,0,0) x-faces
74 // and requires DVom on [-2:nx+1, 0:ny-1] (0,1,0) y-faces
75
76 //
77 // Define UFx, the x-fluxes at cell centers for updating u
78 // (Note that grow arguments are (bx, dir, ng)
79 //
80 ParallelFor(growLo(xbx,0,1),
81 [=] AMREX_GPU_DEVICE (int i, int j, int )
82 {
83 Real uxx_i = ubar(i-1,j,0,krhs)-two*ubar(i ,j,0,krhs)+ubar(i+1,j,0,krhs);
84 Real uxx_ip1 = ubar(i ,j,0,krhs)-two*ubar(i+1,j,0,krhs)+ubar(i+2,j,0,krhs);
85
86 Real Huxx_i = DUon(i-1,j,0)-two*DUon(i ,j,0)+DUon(i+1,j,0);
87 Real Huxx_ip1 = DUon(i ,j,0)-two*DUon(i+1,j,0)+DUon(i+2,j,0);
88
89 if (i == dlo.x && !is_periodic_in_x) {
90 uxx_i = uxx_ip1;
91 Huxx_i = Huxx_ip1;
92 }
93 else if (i == dhi.x && !is_periodic_in_x) {
94 uxx_ip1 = uxx_i;
95 Huxx_ip1 = Huxx_i;
96 }
97
98 Real cff=one/Real(6.0);
99 Real ubar_avg = ubar(i ,j,0,krhs)+ubar(i+1,j,0,krhs);
100
101 UFx(i,j,0)=Real(0.25)*(ubar_avg-cff*(uxx_i+uxx_ip1)) * (DUon(i,j,0)+ DUon(i+1,j,0)-cff*(Huxx_i+ Huxx_ip1));
102 });
103
104 //
105 // Define UFe, the y-fluxes at nodes for updating u
106 //
107 ParallelFor(growHi(xbx,1,1),
108 [=] AMREX_GPU_DEVICE (int i, int j, int )
109 {
110 //should not include grow cells
111 Real uee_j = ubar(i,j-1,0,krhs)-two*ubar(i,j ,0,krhs)+ubar(i,j+1,0,krhs);
112 Real uee_jm1 = ubar(i,j-2,0,krhs)-two*ubar(i,j-1,0,krhs)+ubar(i,j ,0,krhs);
113
114 Real Hvxx_i = DVom(i-1,j,0)-two*DVom(i ,j,0)+DVom(i+1,j,0);
115 Real Hvxx_im1 = DVom(i-2,j,0)-two*DVom(i-1,j,0)+DVom(i ,j,0);
116
117 if (j == dlo.y and !is_periodic_in_y) {
118 uee_jm1 = uee_j;
119 } else if (j == dhi.y+1 and !is_periodic_in_y) {
120 uee_j = uee_jm1;
121 }
122
123 Real cff=one/Real(6.0);
124 Real cff1=ubar(i,j ,0,krhs)+ubar(i,j-1,0,krhs);
125 Real cff2=DVom(i,j,0)+DVom(i-1,j,0);
126
127 UFe(i,j,0)=Real(0.25)*(cff1-(uee_j+uee_jm1)*cff)*
128 (cff2-cff*(Hvxx_i+Hvxx_im1));
129 });
130 } else if (uv_hadv_scheme == AdvectionScheme::centered2) {
131 ParallelFor(growLo(xbx,0,1),
132 [=] AMREX_GPU_DEVICE (int i, int j, int )
133 {
134 UFx(i,j,0) = Real(0.25) * (DUon(i,j,0)+DUon(i+1,j,0)) * (ubar(i,j,0,krhs)+ubar(i+1,j,0,krhs));
135 });
136
137 ParallelFor(growHi(xbx,1,1),
138 [=] AMREX_GPU_DEVICE (int i, int j, int )
139 {
140 UFe(i,j,0) = Real(0.25) * (DVom(i,j,0)+DVom(i-1,j,0)) * (ubar(i,j,0,krhs)+ubar(i,j-1,0,krhs));
141 });
142 }
143
144 //
145 // Add in horizontal advection.
146 //
147 ParallelFor(makeSlab(xbx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
148 {
149 Real cff1=UFx(i,j ,0)-UFx(i-1,j,0);
150 Real cff2=UFe(i,j+1,0)-UFe(i ,j,0);
151
152 rhs_ubar(i,j,0) -= (cff1 + cff2);
153 });
154
155 if (uv_hadv_scheme == AdvectionScheme::upstream3) {
156 // *************************************************************
157 // UPDATING V
158 // *************************************************************
159
160 // Think of the cell-centered box as [ 0:nx-1, 0:ny-1] (0,0,0) cc
161 //
162 // ybx is the y-face-centered box on which we update vbar (with rhs_vbar) [ 0:nx-1, 0:ny ] (1,0,0) y-faces
163 // to do so requires VFe on [ 0:nx-1,-1:ny ] (1,1,0) xy-nodes
164 // which requires vbar on [ 0:nx-1,-2:ny+2] (1,0,0) y-faces
165 // and requires DVom on [ 0:nx-1,-2:ny+2] (0,1,0) x-faces
166 // to do so requires VFx on [ 0:nx , 0:ny ] (0,0,0) cc
167 // which requires vbar on [-2:nx+1, 0:ny ] (1,0,0) y-faces
168 // and requires DUon on [ 0:nx-1,-2:ny+1] (0,0,0) y-faces
169
170 // Grow ybx by one in high x-direction
171 ParallelFor(growHi(ybx,0,1),
172 [=] AMREX_GPU_DEVICE (int i, int j, int )
173 {
174 Real cff=one/Real(6.0);
175 Real vxx_i = vbar(i-1,j,0,krhs)-two*vbar(i ,j,0,krhs)+vbar(i+1,j,0,krhs);
176 Real vxx_im1 = vbar(i-2,j,0,krhs)-two*vbar(i-1,j,0,krhs)+vbar(i ,j,0,krhs);
177
178 Real Huee_j = DUon(i,j-1,0)-two*DUon(i,j ,0)+DUon(i,j+1,0);
179 Real Huee_jm1 = DUon(i,j-2,0)-two*DUon(i,j-1,0)+DUon(i,j ,0);
180
181 if (i == dlo.x and !is_periodic_in_x) {
182 vxx_im1 = vxx_i;
183 } else if (i == dhi.x + 1 and !is_periodic_in_x) {
184 vxx_i = vxx_im1;
185 }
186
187 Real cff1=vbar(i ,j,0,krhs)+vbar(i-1,j,0,krhs);
188 Real cff2=DUon(i,j,0)+DUon(i,j-1,0);
189
190 VFx(i,j,0)=Real(0.25)*(cff1-(vxx_i + vxx_im1)*cff)* (cff2-cff*(Huee_j+ Huee_jm1));
191 });
192
193 ParallelFor(growLo(ybx,1,1),
194 [=] AMREX_GPU_DEVICE (int i, int j, int)
195 {
196 Real vee_j = vbar(i,j-1,0,krhs)-two*vbar(i,j ,0,krhs)+vbar(i,j+1,0,krhs);
197 Real vee_jp1 = vbar(i,j ,0,krhs)-two*vbar(i,j+1,0,krhs)+vbar(i,j+2,0,krhs);
198
199 Real Hvee_j = DVom(i,j-1,0)-two*DVom(i,j ,0)+DVom(i,j+1,0);
200 Real Hvee_jp1 = DVom(i,j ,0)-two*DVom(i,j+1,0)+DVom(i,j+2,0);
201
202 if (j == dlo.y and !is_periodic_in_y) {
203 vee_j = vee_jp1;
204 Hvee_j = Hvee_jp1;
205 }
206 else if (j == dhi.y and !is_periodic_in_y) {
207 vee_jp1 = vee_j;
208 Hvee_jp1 = Hvee_j;
209 }
210
211 Real cff=one/Real(6.0);
212 Real cff1=vbar(i,j ,0,krhs)+vbar(i,j+1,0,krhs);
213
214 VFe(i,j,0) = Real(0.25) * (cff1-(vee_j + vee_jp1)*cff) * (DVom(i,j ,0)+ DVom(i,j+1,0) -
215 cff * (Hvee_j+ Hvee_jp1));
216 });
217 } else if (uv_hadv_scheme == AdvectionScheme::centered2) {
218 ParallelFor(growHi(ybx,0,1),
219 [=] AMREX_GPU_DEVICE (int i, int j, int )
220 {
221 VFx(i,j,0) = Real(0.25) * (DUon(i,j,0) + DUon(i,j-1,0))*(vbar(i,j,0,krhs)+vbar(i-1,j,0,krhs));
222 });
223 ParallelFor(growLo(ybx,1,1),
224 [=] AMREX_GPU_DEVICE (int i, int j, int )
225 {
226 VFe(i,j,0) = Real(0.25) * (DVom(i,j,0) + DVom(i,j+1,0))*(vbar(i,j,0,krhs)+vbar(i,j+1,0,krhs));
227 });
228 }
229
230 ParallelFor(makeSlab(ybx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
231 {
232 Real cff1=VFx(i+1,j,0)-VFx(i ,j,0);
233 Real cff2=VFe(i ,j,0)-VFe(i,j-1,0);
234
235 rhs_vbar(i,j,0) -= (cff1 + cff2);
236 });
237
238}
two
constexpr amrex::Real two
Definition REMORA_Constants.H:8
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
AdvectionScheme::centered2
@ centered2
AdvectionScheme::upstream3
@ upstream3
setVal
mf_h setVal(geomdata.ProbHi(2))
cff1
Real cff1
Definition REMORA_InitAnalyticBiology_BioToy.H:2
cff2
Real cff2
Definition REMORA_InitAnalyticBiology_BioToy.H:3
geomdata
auto geomdata
Definition REMORA_InitAnalyticGridScale_BoundaryLayer.H:2
REMORA::xvel_bc
int xvel_bc() const noexcept
Definition REMORA.H:1313
REMORA::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1578
REMORA::yvel_bc
int yvel_bc() const noexcept
Definition REMORA.H:1314
REMORA::rhs_uv_2d
void rhs_uv_2d(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 const > &Duon, const amrex::Array4< amrex::Real const > &Dvom, const int nrhs)
RHS terms for 2D momentum.
Definition REMORA_rhs_uv_2d.cpp:18
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1717
SolverChoice::uv_Hadv_scheme
AdvectionScheme uv_Hadv_scheme
Definition REMORA_DataStruct.H:909