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