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>(zero);
88 fab_curv.template setVal<RunOn::Device>(zero);
89
90 fab_FX.template setVal<RunOn::Device>(zero);
91 fab_FE.template setVal<RunOn::Device>(zero);
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=one/Real(6.0);
115 Real cffb=one/Real(3.0);
116
117 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
118
119 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
120 {
121 //Upstream3
122 curv(i,j,k)=-FX(i,j,k)+FX(i+1,j,k);
123 });
124
125 //HACK to avoid using the wrong index of t (using upstream3)
126 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
127 {
128 Real max_Huon = std::max(Huon(i,j,k),Real(0.0));
129 Real min_Huon = std::min(Huon(i,j,k),Real(0.0));
130 FX(i,j,k)=Huon(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i-1,j,k))-
131 cffa*(curv(i,j,k)*min_Huon+ curv(i-1,j,k)*max_Huon);
132 });
133
134 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
135
136 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
137 {
138 //Centered4
139 grad(i,j,k)=Real(0.5)*(FX(i,j,k)+FX(i+1,j,k));
140 });
141
142 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
143 {
144 FX(i,j,k)=Huon(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i-1,j,k)-
145 cffb*(grad(i,j,k)- grad(i-1,j,k)));
146 });
147
148 } else {
149 Error("Not a valid horizontal advection scheme");
150 }
151
152 ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
153 {
154 FE(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i,j-1,k,nrhs)) * mskv(i,j,0);
155 });
156
157 Box vtbxp1_slab_lo = makeSlab(vtbxp1,1,dlo.y-1) & vtbxp1;
158 Box vtbxp1_slab_hi = makeSlab(vtbxp1,1,dhi.y+1) & vtbxp1;
159 if (vtbxp1_slab_lo.ok() && !is_periodic_in_y) {
160 ParallelFor(vtbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
161 {
162 FE(i,j,k) = FE(i,j+1,k);
163 });
164 }
165 if (vtbxp1_slab_hi.ok() && !is_periodic_in_y) {
166 ParallelFor(vtbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
167 {
168 FE(i,j+1,k) = FE(i,j,k);
169 });
170 }
171
172
173 cffa=one/Real(6.0);
174 cffb=one/Real(3.0);
175
176 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
177 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
178 {
179 curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k);
180 });
181
182
183 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
184 {
185 Real max_Hvom = std::max(Hvom(i,j,k),Real(0.0));
186 Real min_Hvom = std::min(Hvom(i,j,k),Real(0.0));
187
188 FE(i,j,k)=Hvom(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i,j-1,k))-
189 cffa*(curv(i,j,k)*min_Hvom+ curv(i,j-1,k)*max_Hvom);
190 });
191
192 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
193
194 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
195 {
196 grad(i,j,k)=Real(0.5)*(FE(i,j,k)+FE(i,j+1,k));
197 });
198
199 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
200 {
201 FE(i,j,k)=Hvom(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i,j-1,k)-
202 cffb*(grad(i,j,k)- grad(i,j-1,k)));
203 });
204
205 } else {
206 Error("Not a valid horizontal advection scheme");
207 }
208
209 bool do_rivers_cons = (river_source.size() > 0);
210 if (solverChoice.do_rivers) {
211 int* river_direction_d = river_direction.data();
212 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
213 {
214 int iriver = river_pos(i,j,0);
215 if (iriver >= 0) {
216 if (river_direction_d[iriver] == 0) {
217 if (do_rivers_cons) {
218 FX(i,j,k) = Huon(i,j,k) * river_source(iriver,0,k);
219 } else if ((mskr(i,j,0)==0) && (mskr(i-1,j,0)==1)) {
220 FX(i,j,k) = Huon(i,j,k) * sstore(i-1,j,k);
221 } else if ((mskr(i,j,0)==1) && (mskr(i-1,j,0)==0)) {
222 FX(i,j,k) = Huon(i,j,k) * sstore(i,j,k);
223 }
224 } else if (river_direction_d[iriver] == 1) {
225 if (do_rivers_cons) {
226 FE(i,j,k) = Hvom(i,j,k) * river_source(iriver,0,k);
227 } else if ((mskr(i,j,0)==0) && (mskr(i,j-1,0)==1)) {
228 FE(i,j,k) = Hvom(i,j,k) * sstore(i,j-1,k);
229 } else if ((mskr(i,j,0)==1) && (mskr(i,j-1,0)==0)) {
230 FE(i,j,k) = Hvom(i,j,k) * sstore(i,j,k);
231 }
232 }
233 }
234 });
235 }
236
237 ParallelFor(bx,
238 [=] AMREX_GPU_DEVICE (int i, int j, int k)
239 {
240 //
241 // Add in horizontal advection.
242 //
243 Real cff = dt_lev*pm(i,j,0)*pn(i,j,0);
244 Real cff1=cff*(FX(i+1,j,k)-FX(i,j,k));
245 Real cff2=cff*(FE(i,j+1,k)-FE(i,j,k));
246 Real cff3=cff1+cff2;
247
248 t(i,j,k,nnew) -= cff3;
249 });
250
251 BL_PROFILE_VAR_STOP(phadv);
252 //-----------------------------------------------------------------------
253 // Time-step vertical advection term.
254 //-----------------------------------------------------------------------
255 //Check which type of differences:
256 //
257 // Fourth-order, central differences vertical advective flux
258 // (Tunits m3/s).
259 //
260
261 BL_PROFILE_VAR("REMORA::rhs_t_3d()::vadv",pvadv);
262 ParallelFor(surroundingNodes(bx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
263 {
264 //-----------------------------------------------------------------------
265 // Add in vertical advection.
266 //-----------------------------------------------------------------------
267
268 Real cff1=Real(0.5);
269 Real cff2=Real(7.0)/Real(12.0);
270 Real cff3=one/Real(12.0);
271
272 if (k>=2 && k<=N-1)
273 {
274 FC(i,j,k)=( cff2*(sstore(i ,j,k-1)+ sstore(i,j,k))
275 -cff3*(sstore(i ,j,k-2)+ sstore(i,j,k+1)) ) * ( W(i,j,k));
276 } else if (k==N+1) {
277 FC(i,j,N+1)=Real(0.0);
278 } else if (k==N) {
279 FC(i,j,N)=( cff2*sstore(i ,j,N-1)+ cff1*sstore(i,j,N )
280 -cff3*sstore(i ,j,N-2) ) * ( W(i ,j,N));
281 } else if (k==1) {
282 FC(i,j,1)=( cff2*sstore(i ,j,1)+ cff1*sstore(i,j,0)
283 -cff3*sstore(i ,j,2) ) * ( W(i ,j,1));
284 } else if (k==0) {
285 FC(i,j,0) = Real(0.0);
286 }
287 });
288
289 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
290 {
291 Real cff1=dt_lev*pm(i,j,0)*pn(i,j,0);
292 Real cff4=FC(i,j,k+1)-FC(i,j,k);
293
294 t(i,j,k) = (t(i,j,k)-cff1*cff4) / Hz(i,j,k);
295 });
296 BL_PROFILE_VAR_STOP(pvadv);
297}
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
AdvectionScheme::upstream3
@ upstream3
AdvectionScheme::centered4
@ centered4
NGROW
#define NGROW
Definition REMORA_IndexDefines.H:30
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::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:1717
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:1496
SolverChoice::tracer_Hadv_scheme
AdvectionScheme tracer_Hadv_scheme
Definition REMORA_DataStruct.H:908