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 const Array4<Real>& fx_reg,
45 const Array4<Real>& fy_reg)
46{
47 BL_PROFILE("REMORA::rhs_t_3d()");
48 const Box& domain = geom[lev].Domain();
49 const auto dlo = amrex::lbound(domain);
50 const auto dhi = amrex::ubound(domain);
51
52 GeometryData const& geomdata = geom[0].data();
53 bool is_periodic_in_x = geomdata.isPeriodic(0);
54 bool is_periodic_in_y = geomdata.isPeriodic(1);
55
56 //copy the tilebox
57 Box tbxp1x = bx;
58 Box tbxp1y = bx;
59 Box tbxp1 = bx;
60 Box tbxp2 = bx;
61
62 tbxp2.grow(IntVect(NGROW,NGROW,0));
63 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
64 tbxp1x.grow(IntVect(NGROW-1,0,0));
65 tbxp1y.grow(IntVect(0,NGROW-1,0));
66
67 // Because grad, curv, FX, FE, are all local, do surroundinNodes
68 Box utbxp1 = surroundingNodes(tbxp1x, 0);
69 Box vtbxp1 = surroundingNodes(tbxp1y, 1);
70 Box ubx = surroundingNodes(bx, 0);
71 Box vbx = surroundingNodes(bx, 1);
72
73
74 //
75 // Scratch space
76 //
77 FArrayBox fab_grad(tbxp2,1,amrex::The_Async_Arena());
78 FArrayBox fab_curv(tbxp2,1,amrex::The_Async_Arena());
79
80 FArrayBox fab_FX(tbxp2,1,amrex::The_Async_Arena());
81 FArrayBox fab_FE(tbxp2,1,amrex::The_Async_Arena());
82
83 auto curv=fab_curv.array();
84 auto grad=fab_grad.array();
85
86 auto FX=fab_FX.array();
87 auto FE=fab_FE.array();
88
89 fab_grad.template setVal<RunOn::Device>(zero);
90 fab_curv.template setVal<RunOn::Device>(zero);
91
92 fab_FX.template setVal<RunOn::Device>(zero);
93 fab_FE.template setVal<RunOn::Device>(zero);
94
95 BL_PROFILE_VAR("REMORA::rhs_t_3d()::hadv",phadv);
96 ParallelFor(utbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
97 {
98 FX(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i-1,j,k,nrhs)) * msku(i,j,0);
99 });
100
101 Box utbxp1_slab_lo = makeSlab(utbxp1,0,dlo.x-1) & utbxp1;
102 Box utbxp1_slab_hi = makeSlab(utbxp1,0,dhi.x+1) & utbxp1;
103 if (utbxp1_slab_lo.ok() && !is_periodic_in_x) {
104 ParallelFor(utbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
105 {
106 FX(i,j,k) = FX(i+1,j,k);
107 });
108 }
109 if (utbxp1_slab_hi.ok() && !is_periodic_in_x) {
110 ParallelFor(utbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
111 {
112 FX(i+1,j,k) = FX(i,j,k);
113 });
114 }
115
116 Real cffa=one/Real(6.0);
117 Real cffb=one/Real(3.0);
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 //HACK to avoid using the wrong index of t (using upstream3)
128 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
129 {
130 Real max_Huon = std::max(Huon(i,j,k),Real(0.0));
131 Real min_Huon = std::min(Huon(i,j,k),Real(0.0));
132 FX(i,j,k)=Huon(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i-1,j,k))-
133 cffa*(curv(i,j,k)*min_Huon+ curv(i-1,j,k)*max_Huon);
134 });
135
136 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
137
138 ParallelFor(tbxp1x, [=] AMREX_GPU_DEVICE (int i, int j, int k)
139 {
140 //Centered4
141 grad(i,j,k)=Real(0.5)*(FX(i,j,k)+FX(i+1,j,k));
142 });
143
144 ParallelFor(ubx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
145 {
146 FX(i,j,k)=Huon(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i-1,j,k)-
147 cffb*(grad(i,j,k)- grad(i-1,j,k)));
148 });
149
150 } else {
151 Error("Not a valid horizontal advection scheme");
152 }
153
154 ParallelFor(vtbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
155 {
156 FE(i,j,k)=(sstore(i,j,k,nrhs)-sstore(i,j-1,k,nrhs)) * mskv(i,j,0);
157 });
158
159 Box vtbxp1_slab_lo = makeSlab(vtbxp1,1,dlo.y-1) & vtbxp1;
160 Box vtbxp1_slab_hi = makeSlab(vtbxp1,1,dhi.y+1) & vtbxp1;
161 if (vtbxp1_slab_lo.ok() && !is_periodic_in_y) {
162 ParallelFor(vtbxp1_slab_lo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
163 {
164 FE(i,j,k) = FE(i,j+1,k);
165 });
166 }
167 if (vtbxp1_slab_hi.ok() && !is_periodic_in_y) {
168 ParallelFor(vtbxp1_slab_hi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
169 {
170 FE(i,j+1,k) = FE(i,j,k);
171 });
172 }
173
174
175 cffa=one/Real(6.0);
176 cffb=one/Real(3.0);
177
178 if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::upstream3) {
179 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
180 {
181 curv(i,j,k)=-FE(i,j,k)+FE(i,j+1,k);
182 });
183
184
185 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
186 {
187 Real max_Hvom = std::max(Hvom(i,j,k),Real(0.0));
188 Real min_Hvom = std::min(Hvom(i,j,k),Real(0.0));
189
190 FE(i,j,k)=Hvom(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i,j-1,k))-
191 cffa*(curv(i,j,k)*min_Hvom+ curv(i,j-1,k)*max_Hvom);
192 });
193
194 } else if (solverChoice.tracer_Hadv_scheme == AdvectionScheme::centered4) {
195
196 ParallelFor(tbxp1y, [=] AMREX_GPU_DEVICE (int i, int j, int k)
197 {
198 grad(i,j,k)=Real(0.5)*(FE(i,j,k)+FE(i,j+1,k));
199 });
200
201 ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
202 {
203 FE(i,j,k)=Hvom(i,j,k)*Real(0.5)*(sstore(i,j,k)+sstore(i,j-1,k)-
204 cffb*(grad(i,j,k)- grad(i,j-1,k)));
205 });
206
207 } else {
208 Error("Not a valid horizontal advection scheme");
209 }
210
211 bool do_rivers_cons = (river_source.size() > 0);
212 if (solverChoice.do_rivers) {
213 int* river_direction_d = river_direction.data();
214 ParallelFor(tbxp1, [=] AMREX_GPU_DEVICE (int i, int j, int k)
215 {
216 int iriver = river_pos(i,j,0);
217 if (iriver >= 0) {
218 if (river_direction_d[iriver] == 0) {
219 if (do_rivers_cons) {
220 FX(i,j,k) = Huon(i,j,k) * river_source(iriver,0,k);
221 } else if ((mskr(i,j,0)==0) && (mskr(i-1,j,0)==1)) {
222 FX(i,j,k) = Huon(i,j,k) * sstore(i-1,j,k);
223 } else if ((mskr(i,j,0)==1) && (mskr(i-1,j,0)==0)) {
224 FX(i,j,k) = Huon(i,j,k) * sstore(i,j,k);
225 }
226 } else if (river_direction_d[iriver] == 1) {
227 if (do_rivers_cons) {
228 FE(i,j,k) = Hvom(i,j,k) * river_source(iriver,0,k);
229 } else if ((mskr(i,j,0)==0) && (mskr(i,j-1,0)==1)) {
230 FE(i,j,k) = Hvom(i,j,k) * sstore(i,j-1,k);
231 } else if ((mskr(i,j,0)==1) && (mskr(i,j-1,0)==0)) {
232 FE(i,j,k) = Hvom(i,j,k) * sstore(i,j,k);
233 }
234 }
235 }
236 });
237 }
238
239 ParallelFor(bx,
240 [=] AMREX_GPU_DEVICE (int i, int j, int k)
241 {
242 //
243 // Add in horizontal advection.
244 //
245 Real cff = dt_lev*pm(i,j,0)*pn(i,j,0);
246 Real cff1=cff*(FX(i+1,j,k)-FX(i,j,k));
247 Real cff2=cff*(FE(i,j+1,k)-FE(i,j,k));
248 Real cff3=cff1+cff2;
249
250 t(i,j,k,nnew) -= cff3;
251 });
252
253 // Hand the same fluxes to the flux register, scaled so its internal dt/dx reproduces
254 // the dt*pm*pn above. REMORA's divergence uses the receiving cell's metric, not a face
255 // one, so the face weight is an average: exact where pm and pn are uniform.
256 if (fx_reg) {
257 const Real dx0 = geom[lev].CellSize(0);
258 const Real dx1 = geom[lev].CellSize(1);
259 ParallelFor(surroundingNodes(bx,0), [=] AMREX_GPU_DEVICE (int i, int j, int k)
260 {
261 Real w = Real(0.5) * (pm(i-1,j,0)*pn(i-1,j,0) + pm(i,j,0)*pn(i,j,0));
262 fx_reg(i,j,k) = FX(i,j,k) * w * dx0;
263 });
264 ParallelFor(surroundingNodes(bx,1), [=] AMREX_GPU_DEVICE (int i, int j, int k)
265 {
266 Real w = Real(0.5) * (pm(i,j-1,0)*pn(i,j-1,0) + pm(i,j,0)*pn(i,j,0));
267 fy_reg(i,j,k) = FE(i,j,k) * w * dx1;
268 });
269 }
270
271 BL_PROFILE_VAR_STOP(phadv);
272 //-----------------------------------------------------------------------
273 // Time-step vertical advection term.
274 //-----------------------------------------------------------------------
275 //Check which type of differences:
276 //
277 // Fourth-order, central differences vertical advective flux
278 // (Tunits m3/s).
279 //
280
281 BL_PROFILE_VAR("REMORA::rhs_t_3d()::vadv",pvadv);
282 ParallelFor(surroundingNodes(bx,2), [=] AMREX_GPU_DEVICE (int i, int j, int k)
283 {
284 //-----------------------------------------------------------------------
285 // Add in vertical advection.
286 //-----------------------------------------------------------------------
287
288 Real cff1=Real(0.5);
289 Real cff2=Real(7.0)/Real(12.0);
290 Real cff3=one/Real(12.0);
291
292 if (k>=2 && k<=N-1)
293 {
294 FC(i,j,k)=( cff2*(sstore(i ,j,k-1)+ sstore(i,j,k))
295 -cff3*(sstore(i ,j,k-2)+ sstore(i,j,k+1)) ) * ( W(i,j,k));
296 } else if (k==N+1) {
297 FC(i,j,N+1)=Real(0.0);
298 } else if (k==N) {
299 FC(i,j,N)=( cff2*sstore(i ,j,N-1)+ cff1*sstore(i,j,N )
300 -cff3*sstore(i ,j,N-2) ) * ( W(i ,j,N));
301 } else if (k==1) {
302 FC(i,j,1)=( cff2*sstore(i ,j,1)+ cff1*sstore(i,j,0)
303 -cff3*sstore(i ,j,2) ) * ( W(i ,j,1));
304 } else if (k==0) {
305 FC(i,j,0) = Real(0.0);
306 }
307 });
308
309 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
310 {
311 Real cff1=dt_lev*pm(i,j,0)*pn(i,j,0);
312 Real cff4=FC(i,j,k+1)-FC(i,j,k);
313
314 t(i,j,k) = (t(i,j,k)-cff1*cff4) / Hz(i,j,k);
315 });
316 BL_PROFILE_VAR_STOP(pvadv);
317}
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_InitAnalyticProb_Advection.H:12
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, const amrex::Array4< amrex::Real > &fx_reg=amrex::Array4< amrex::Real >(), const amrex::Array4< amrex::Real > &fy_reg=amrex::Array4< amrex::Real >())
RHS terms for tracer.
Definition REMORA_rhs_t_3d.cpp:28
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1949
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:1629
SolverChoice::tracer_Hadv_scheme
AdvectionScheme tracer_Hadv_scheme
Definition REMORA_DataStruct.H:924