REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_rhs_t_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 ] bx tilebox
8 * @param[inout] t tracer data
9 * @param[in ] sstore scratch space for tracer calculations
10 * @param[in ] Huon u-volume flux
11 * @param[in ] Hvom v-volume flux
12 * @param[in ] Hz vertical cell height
13 * @param[in ] pn 1/dx
14 * @param[in ] pm 1/dy
15 * @param[in ] W vertical velocity
16 * @param FC temporary
17 * @param[in ] mskr land-sea mask on rho-points
18 * @param[in ] msku land-sea mask on u-points
19 * @param[in ] mskv land-sea mask on v-points
20 * @param[in ] river_pos river positions
21 * @param[in ] river_source river source data to add, if using
22 * @param[in ] nrhs index of RHS component
23 * @param[in ] nnew index of current time step
24 * @param[in ] N number of vertical levels
25 * @param[in ] dt_lev time step at this level
26 */
27void
28REMORA::rhs_t_3d (int lev, const Box& bx,
29 const Array4<Real >& t,
30 const Array4<Real const>& sstore,
31 const Array4<Real const>& Huon,
32 const Array4<Real const>& Hvom,
33 const Array4<Real const>& Hz,
34 const Array4<Real const>& pn,
35 const Array4<Real const>& pm,
36 const Array4<Real const>& W ,
37 const Array4<Real >& FC,
38 const Array4<Real const>& mskr,
39 const Array4<Real const>& msku,
40 const Array4<Real const>& mskv,
41 const Array4<int const>& river_pos,
42 const Array4<Real const>& river_source,
43 int nrhs, int nnew, int N, Real dt_lev)
44{
45 BL_PROFILE("REMORA::rhs_t_3d()");
46 const Box& domain = geom[lev].Domain();
47 const auto dlo = amrex::lbound(domain);
48 const auto dhi = amrex::ubound(domain);
49
50 GeometryData const& geomdata = geom[0].data();
51 bool is_periodic_in_x = geomdata.isPeriodic(0);
52 bool is_periodic_in_y = geomdata.isPeriodic(1);
53
54 //copy the tilebox
55 Box tbxp1x = bx;
56 Box tbxp1y = bx;
57 Box tbxp1 = bx;
58 Box tbxp2 = bx;
59
60 tbxp2.grow(IntVect(NGROW,NGROW,0));
61 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
62 tbxp1x.grow(IntVect(NGROW-1,0,0));
63 tbxp1y.grow(IntVect(0,NGROW-1,0));
64
65 // Because grad, curv, FX, FE, are all local, do surroundinNodes
66 Box utbxp1 = surroundingNodes(tbxp1x, 0);
67 Box vtbxp1 = surroundingNodes(tbxp1y, 1);
68 Box ubx = surroundingNodes(bx, 0);
69 Box vbx = surroundingNodes(bx, 1);
70
71
72 //
73 // Scratch space
74 //
75 FArrayBox fab_grad(tbxp2,1,amrex::The_Async_Arena());
76 FArrayBox fab_curv(tbxp2,1,amrex::The_Async_Arena());
77
78 FArrayBox fab_FX(tbxp2,1,amrex::The_Async_Arena());
79 FArrayBox fab_FE(tbxp2,1,amrex::The_Async_Arena());
80
81 auto curv=fab_curv.array();
82 auto grad=fab_grad.array();
83
84 auto FX=fab_FX.array();
85 auto FE=fab_FE.array();
86
87 fab_grad.template setVal<RunOn::Device>(0.);
88 fab_curv.template setVal<RunOn::Device>(0.);
89
90 fab_FX.template setVal<RunOn::Device>(0.);
91 fab_FE.template setVal<RunOn::Device>(0.);
92
93 BL_PROFILE_VAR("REMORA::rhs_t_3d()::hadv",phadv);
94 ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
95 {
96 FX(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i-1,j,k,nrhs)) * msku(i,j,0);
97 });
98
99 Box utbxp1_slab_lo = makeSlab(utbxp1,0,dlo.x-1) & utbxp1;
100 Box utbxp1_slab_hi = makeSlab(utbxp1,0,dhi.x+1) & utbxp1;
101 if (utbxp1_slab_lo.ok() && !is_periodic_in_x) {
102 ParallelFor(utbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
103 {
104 FX(i,j,k) = FX(i+1,j,k);
105 });
106 }
107 if (utbxp1_slab_hi.ok() && !is_periodic_in_x) {
108 ParallelFor(utbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
109 {
110 FX(i+1,j,k) = FX(i,j,k);
111 });
112 }
113
114 Real cffa=1.0_rt/6.0_rt;
115 Real cffb=1.0_rt/3.0_rt;
116
117 if(solverChoice.flat_bathymetry) {
118
119 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
120
121 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
122 {
123 //Upstream3
124 curv(i,j,k)=-FX(i,j,k)+FX(i+1,j,k);
125 });
126
127 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
128 {
129 Real max_Huon = std::max(Huon(i,j,k),0.0_rt);
130 Real min_Huon = std::min(Huon(i,j,k),0.0_rt);
131 FX(i,j,k)=Huon(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i-1,j,k))+
132 cffa*(curv(i,j,k)*min_Huon+ curv(i-1,j,k)*max_Huon);
133 });
134
135 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
136
137 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
138 {
139 //Centered4
140 grad(i,j,k)=0.5_rt*(FX(i,j,k)+FX(i+1,j,k));
141 });
142
143 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
144 {
145 FX(i,j,k)=Huon(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i-1,j,k))+
146 cffb*(grad(i,j,k)+ grad(i-1,j,k));
147 });
148
149 } else {
150 Error("Not a valid horizontal advection scheme");
151 }
152
153 } else {
154
155 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
156
157 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
158 {
159 //Upstream3
160 curv(i,j,k)=-FX(i,j,k)+FX(i+1,j,k);
161 });
162
163 //HACK to avoid using the wrong index of t (using upstream3)
164 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
165 {
166 Real max_Huon = std::max(Huon(i,j,k),0.0_rt);
167 Real min_Huon = std::min(Huon(i,j,k),0.0_rt);
168 FX(i,j,k)=Huon(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i-1,j,k))-
169 cffa*(curv(i,j,k)*min_Huon+ curv(i-1,j,k)*max_Huon);
170 });
171
172 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
173
174 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
175 {
176 //Centered4
177 grad(i,j,k)=0.5_rt*(FX(i,j,k)+FX(i+1,j,k));
178 });
179
180 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
181 {
182 FX(i,j,k)=Huon(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i-1,j,k)-
183 cffb*(grad(i,j,k)- grad(i-1,j,k)));
184 });
185
186 } else {
187 Error("Not a valid horizontal advection scheme");
188 }
189 } // flat bathymetry?
190
191 ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
192 {
193 FE(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i,j-1,k,nrhs)) * mskv(i,j,0);
194 });
195
196 Box vtbxp1_slab_lo = makeSlab(vtbxp1,1,dlo.y-1) & vtbxp1;
197 Box vtbxp1_slab_hi = makeSlab(vtbxp1,1,dhi.y+1) & vtbxp1;
198 if (vtbxp1_slab_lo.ok() && !is_periodic_in_y) {
199 ParallelFor(vtbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
200 {
201 FE(i,j,k) = FE(i,j+1,k);
202 });
203 }
204 if (vtbxp1_slab_hi.ok() && !is_periodic_in_y) {
205 ParallelFor(vtbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
206 {
207 FE(i,j+1,k) = FE(i,j,k);
208 });
209 }
210
211
212 cffa=1.0_rt/6.0_rt;
213 cffb=1.0_rt/3.0_rt;
214 if (solverChoice.flat_bathymetry) {
215
216 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
217
218 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
219 {
220 curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k);
221 });
222
223 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
224 {
225 Real max_Hvom = std::max(Hvom(i,j,k),0.0_rt);
226 Real min_Hvom = std::min(Hvom(i,j,k),0.0_rt);
227
228 FE(i,j,k)=Hvom(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i,j-1,k))+
229 cffa*(curv(i,j,k)*min_Hvom+ curv(i,j-1,k)*max_Hvom);
230 });
231
232 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
233
234 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
235 {
236 grad(i,j,k)=0.5_rt*(FE(i,j,k)+FE(i,j+1,k));
237 });
238
239 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
240 {
241 FE(i,j,k)=Hvom(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i,j-1,k))+
242 cffb*(grad(i,j,k)+ grad(i,j-1,k));
243 });
244
245 } else {
246 Error("Not a valid horizontal advection scheme");
247 }
248
249 } else {
250
251 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
252 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
253 {
254 curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k);
255 });
256
257
258 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
259 {
260 Real max_Hvom = std::max(Hvom(i,j,k),0.0_rt);
261 Real min_Hvom = std::min(Hvom(i,j,k),0.0_rt);
262
263 FE(i,j,k)=Hvom(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i,j-1,k))-
264 cffa*(curv(i,j,k)*min_Hvom+ curv(i,j-1,k)*max_Hvom);
265 });
266
267 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
268
269 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
270 {
271 grad(i,j,k)=0.5_rt*(FE(i,j,k)+FE(i,j+1,k));
272 });
273
274 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
275 {
276 FE(i,j,k)=Hvom(i,j,k)*0.5_rt*(sstore(i,j,k)+sstore(i,j-1,k)-
277 cffb*(grad(i,j,k)- grad(i,j-1,k)));
278 });
279
280 } else {
281 Error("Not a valid horizontal advection scheme");
282 }
283 }
284
285 bool do_rivers_cons = (river_source.size() > 0);
286 if (solverChoice.do_rivers) {
287 int* river_direction_d = river_direction.data();
288 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
289 {
290 int iriver = river_pos(i,j,0);
291 if (iriver >= 0) {
292 if (river_direction_d[iriver] == 0) {
293 if (do_rivers_cons) {
294 FX(i,j,k) = Huon(i,j,k) * river_source(iriver,0,k);
295 } else if ((mskr(i,j,0)==0) && (mskr(i-1,j,0)==1)) {
296 FX(i,j,k) = Huon(i,j,k) * sstore(i-1,j,k);
297 } else if ((mskr(i,j,0)==1) && (mskr(i-1,j,0)==0)) {
298 FX(i,j,k) = Huon(i,j,k) * sstore(i,j,k);
299 }
300 } else if (river_direction_d[iriver] == 1) {
301 if (do_rivers_cons) {
302 FE(i,j,k) = Hvom(i,j,k) * river_source(iriver,0,k);
303 } else if ((mskr(i,j,0)==0) && (mskr(i,j-1,0)==1)) {
304 FE(i,j,k) = Hvom(i,j,k) * sstore(i,j-1,k);
305 } else if ((mskr(i,j,0)==1) && (mskr(i,j-1,0)==0)) {
306 FE(i,j,k) = Hvom(i,j,k) * sstore(i,j,k);
307 }
308 }
309 }
310 });
311 }
312
313 ParallelFor(bx,
314 [=] AMREX_GPU_DEVICE (int i, int j, int k)
315 {
316 //
317 // Add in horizontal advection.
318 //
319 Real cff = dt_lev*pm(i,j,0)*pn(i,j,0);
320 Real cff1=cff*(FX(i+1,j,k)-FX(i,j,k));
321 Real cff2=cff*(FE(i,j+1,k)-FE(i,j,k));
322 Real cff3=cff1+cff2;
323
324 t(i,j,k,nnew) -= cff3;
325 });
326
327 BL_PROFILE_VAR_STOP(phadv);
328 //-----------------------------------------------------------------------
329 // Time-step vertical advection term.
330 //-----------------------------------------------------------------------
331 //Check which type of differences:
332 //
333 // Fourth-order, central differences vertical advective flux
334 // (Tunits m3/s).
335 //
336
337 BL_PROFILE_VAR("REMORA::rhs_t_3d()::vadv",pvadv);
338 ParallelFor(surroundingNodes(bx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
339 {
340 //-----------------------------------------------------------------------
341 // Add in vertical advection.
342 //-----------------------------------------------------------------------
343
344 Real cff1=0.5_rt;
345 Real cff2=7.0_rt/12.0_rt;
346 Real cff3=1.0_rt/12.0_rt;
347
348 if (k>=2 && k<=N-1)
349 {
350 FC(i,j,k)=( cff2*(sstore(i ,j,k-1)+ sstore(i,j,k))
351 -cff3*(sstore(i ,j,k-2)+ sstore(i,j,k+1)) ) * ( W(i,j,k));
352 } else if (k==N+1) {
353 FC(i,j,N+1)=0.0_rt;
354 } else if (k==N) {
355 FC(i,j,N)=( cff2*sstore(i ,j,N-1)+ cff1*sstore(i,j,N )
356 -cff3*sstore(i ,j,N-2) ) * ( W(i ,j,N));
357 } else if (k==1) {
358 FC(i,j,1)=( cff2*sstore(i ,j,1)+ cff1*sstore(i,j,0)
359 -cff3*sstore(i ,j,2) ) * ( W(i ,j,1));
360 } else if (k==0) {
361 FC(i,j,0) = 0.0_rt;
362 }
363 });
364
365 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
366 {
367 Real cff1=dt_lev*pm(i,j,0)*pn(i,j,0);
368 Real cff4=FC(i,j,k+1)-FC(i,j,k);
369
370 t(i,j,k) = (t(i,j,k)-cff1*cff4) / Hz(i,j,k);
371 });
372 BL_PROFILE_VAR_STOP(pvadv);
373}
AdvectionScheme::upstream3
@ upstream3
AdvectionScheme::centered4
@ centered4
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:13
REMORA::rhs_t_3d
void rhs_t_3d(int lev, const amrex::Box &bx, const amrex::Array4< amrex::Real > &t, const amrex::Array4< amrex::Real const > &tempstore, const amrex::Array4< amrex::Real const > &Huon, const amrex::Array4< amrex::Real const > &Hvom, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &W, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &mskr, const amrex::Array4< amrex::Real const > &msku, const amrex::Array4< amrex::Real const > &mskv, const amrex::Array4< int const > &river_pos, const amrex::Array4< amrex::Real const > &river_source, int nrhs, int nnew, int N, const amrex::Real dt_lev)
RHS terms for tracer.
Definition REMORA_rhs_t_3d.cpp:28
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
REMORA::river_direction
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1096
SolverChoice::tracer_Hadv_scheme
AdvectionScheme tracer_Hadv_scheme
Definition REMORA_DataStruct.H:571