AMReX-Hydro
AMReX-based hydro routines for low Mach number flows
 
Loading...
Searching...
No Matches
hydro_ebmol_edge_state_K.H
Go to the documentation of this file.
1/**
2 * \file hydro_ebmol_edge_state_K.H
3 * \addtogroup EBMOL
4 * @{
5 *
6 */
7
8#ifndef HYDRO_EBMOL_EDGE_STATE_K_H_
9#define HYDRO_EBMOL_EDGE_STATE_K_H_
10
11#include <AMReX_EB_Slopes_K.H>
12#include <AMReX_BCRec.H>
13#include <hydro_constants.H>
14#include <hydro_bcs_K.H>
15
16namespace EBMOL {
17
18AMREX_GPU_DEVICE AMREX_FORCE_INLINE
19amrex::Real hydro_ebmol_xedge_state_extdir ( AMREX_D_DECL(int i, int j, int k), int n,
20 amrex::Array4<amrex::Real const> const& q,
21 amrex::Array4<amrex::Real const> const& umac,
22 AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
23 amrex::Array4<amrex::Real const> const& fcy,
24 amrex::Array4<amrex::Real const> const& fcz),
25 amrex::Array4<amrex::Real const> const& ccc,
26 amrex::Array4<amrex::Real const> const& vfrac,
27 amrex::Array4<amrex::EBCellFlag const> const& flag,
28 amrex::BCRec const* const d_bcrec,
29 amrex::Box const& domain,
30 int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
31{
32
33#if (AMREX_SPACEDIM==2)
34 const int k = 0;
35#endif
36
37 amrex::Real qs;
38 int domlo = domain.smallEnd(0);
39 int domhi = domain.bigEnd(0);
40 bool extlo = d_bcrec[n].lo(0) == amrex::BCType::ext_dir;
41 bool exthi = d_bcrec[n].hi(0) == amrex::BCType::ext_dir;
42
43 // Only short-circuit on faces outside the domain. On the domain face itself we
44 // fall through to SetEdgeBCsLo/Hi, which keeps the Dirichlet value for the
45 // inflowing state but lets the interior state be transported out, so that an
46 // ext_dir face with local backflow upwinds the same way (EB)Godunov does.
47 if ( extlo && i < domlo)
48 {
49 qs = q(domlo-1,j,k,n);
50 }
51 else if ( exthi && i > domhi+1)
52 {
53 qs = q(domhi+1,j,k,n);
54 }
55 else
56 {
57 const int domain_ilo = domain.smallEnd(0);
58 const int domain_ihi = domain.bigEnd(0);
59 const int domain_jlo = domain.smallEnd(1);
60 const int domain_jhi = domain.bigEnd(1);
61#if (AMREX_SPACEDIM == 3)
62 const int domain_klo = domain.smallEnd(2);
63 const int domain_khi = domain.bigEnd(2);
64#endif
65
66 amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of z-face we are extrapolating to
67#if (AMREX_SPACEDIM==3)
68 amrex::Real zf = fcx(i,j,k,1);
69#endif
70
71 AMREX_D_TERM(amrex::Real xc = ccc(i,j,k,0);, // centroid of cell (i,j,k)
72 amrex::Real yc = ccc(i,j,k,1);,
73 amrex::Real zc = ccc(i,j,k,2););
74
75 AMREX_D_TERM(amrex::Real delta_x = amrex::Real(0.5) + xc;,
76 amrex::Real delta_y = yf - yc;,
77 amrex::Real delta_z = zf - zc;);
78
79 amrex::Real cc_qmax = amrex::max(q(i,j,k,n),q(i-1,j,k,n));
80 amrex::Real cc_qmin = amrex::min(q(i,j,k,n),q(i-1,j,k,n));
81
82 AMREX_D_TERM(bool extdir_or_ho_ilo = (d_bcrec[n].lo(0) == amrex::BCType::ext_dir) ||
83 (d_bcrec[n].lo(0) == amrex::BCType::hoextrap);,
84 bool extdir_or_ho_jlo = (d_bcrec[n].lo(1) == amrex::BCType::ext_dir) ||
85 (d_bcrec[n].lo(1) == amrex::BCType::hoextrap);,
86 bool extdir_or_ho_klo = (d_bcrec[n].lo(2) == amrex::BCType::ext_dir) ||
87 (d_bcrec[n].lo(2) == amrex::BCType::hoextrap););
88
89 AMREX_D_TERM(bool extdir_or_ho_ihi = (d_bcrec[n].hi(0) == amrex::BCType::ext_dir) ||
90 (d_bcrec[n].hi(0) == amrex::BCType::hoextrap);,
91 bool extdir_or_ho_jhi = (d_bcrec[n].hi(1) == amrex::BCType::ext_dir) ||
92 (d_bcrec[n].hi(1) == amrex::BCType::hoextrap);,
93 bool extdir_or_ho_khi = (d_bcrec[n].hi(2) == amrex::BCType::ext_dir) ||
94 (d_bcrec[n].hi(2) == amrex::BCType::hoextrap););
95
96 const auto& slopes_eb_hi =
97 amrex_lim_slopes_extdir_eb(i, j, k, n, q, ccc, vfrac,
98 AMREX_D_DECL(fcx,fcy,fcz), flag,
99 AMREX_D_DECL(extdir_or_ho_ilo, extdir_or_ho_jlo, extdir_or_ho_klo),
100 AMREX_D_DECL(extdir_or_ho_ihi, extdir_or_ho_jhi, extdir_or_ho_khi),
101 AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
102 AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
103 order);
104
105
106#if (AMREX_SPACEDIM==3)
107 amrex::Real qpls = q(i ,j,k,n) - delta_x * slopes_eb_hi[0]
108 + delta_y * slopes_eb_hi[1]
109 + delta_z * slopes_eb_hi[2];
110#else
111 amrex::Real qpls = q(i ,j,k,n) - delta_x * slopes_eb_hi[0]
112 + delta_y * slopes_eb_hi[1];
113#endif
114
115 qpls = amrex::max(amrex::min(qpls, cc_qmax), cc_qmin);
116
117 AMREX_D_TERM(xc = ccc(i-1,j,k,0);, // centroid of cell (i-1,j,k)
118 yc = ccc(i-1,j,k,1);,
119 zc = ccc(i-1,j,k,2););
120
121 AMREX_D_TERM( delta_x = amrex::Real(0.5) - xc;,
122 delta_y = yf - yc;,
123 delta_z = zf - zc;);
124
125 // Compute slopes of component "n" of q
126 const auto& slopes_eb_lo =
127 amrex_lim_slopes_extdir_eb(i-1, j, k, n, q, ccc, vfrac,
128 AMREX_D_DECL(fcx,fcy,fcz), flag,
129 AMREX_D_DECL(extdir_or_ho_ilo, extdir_or_ho_jlo, extdir_or_ho_klo),
130 AMREX_D_DECL(extdir_or_ho_ihi, extdir_or_ho_jhi, extdir_or_ho_khi),
131 AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
132 AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
133 order);
134
135#if (AMREX_SPACEDIM==3)
136 amrex::Real qmns = q(i-1,j,k,n) + delta_x * slopes_eb_lo[0]
137 + delta_y * slopes_eb_lo[1]
138 + delta_z * slopes_eb_lo[2];
139#else
140 amrex::Real qmns = q(i-1,j,k,n) + delta_x * slopes_eb_lo[0]
141 + delta_y * slopes_eb_lo[1];
142#endif
143
144 qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
145
146 HydroBC::SetEdgeBCsLo(0, i, j, k, n, q, qmns, qpls, umac(i,j,k), d_bcrec[n].lo(0), domain_ilo, is_velocity);
147 HydroBC::SetEdgeBCsHi(0, i, j, k, n, q, qmns, qpls, umac(i,j,k), d_bcrec[n].hi(0), domain_ihi, is_velocity);
148
149 // The clamp below forces the normal velocity to be outflowing at an outflow
150 // face. allow_inflow_on_outflow lets the caller keep an inflowing MAC
151 // velocity there, so the state must not be clamped against it.
152 if (!allow_inflow_on_outflow)
153 {
154 if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
155 {
156 if ( umac(i,j,k) >= amrex::Real(0.) && n==XVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0.0));
157 qmns = qpls;
158 }
159 if ( (i==domain_ihi+1) && (d_bcrec[n].hi(0) == amrex::BCType::foextrap || d_bcrec[n].hi(0) == amrex::BCType::hoextrap) )
160 {
161 if ( umac(i,j,k) <= amrex::Real(0.) && n==XVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0.0));
162 qpls = qmns;
163 }
164 }
165
166 if (umac(i,j,k) > small_vel)
167 {
168 qs = qmns;
169 }
170 else if (umac(i,j,k) < - small_vel)
171 {
172 qs = qpls;
173 }
174 else
175 {
176 qs = amrex::Real(0.5)*(qmns+qpls);
177 }
178 }
179
180 return qs;
181}
182
183AMREX_GPU_DEVICE AMREX_FORCE_INLINE
184amrex::Real hydro_ebmol_xedge_state ( AMREX_D_DECL(int i, int j, int k), int n,
185 amrex::Array4<amrex::Real const> const& q,
186 amrex::Array4<amrex::Real const> const& umac,
187 AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
188 amrex::Array4<amrex::Real const> const& fcy,
189 amrex::Array4<amrex::Real const> const& fcz),
190 amrex::Array4<amrex::Real const> const& ccc,
191 amrex::Array4<amrex::Real const> const& vfrac,
192 amrex::Array4<amrex::EBCellFlag const> const& flag,
193 amrex::BCRec const* const d_bcrec,
194 amrex::Box const& domain,
195 int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
196{
197#if (AMREX_SPACEDIM==2)
198 const int k = 0;
199#endif
200
201 const int domain_ilo = domain.smallEnd(0);
202 const int domain_ihi = domain.bigEnd(0);
203
204 // local (y,z) of centroid of x-face we are extrapolating to
205 amrex::Real yf = fcx(i,j,k,0);
206#if (AMREX_SPACEDIM==3)
207 amrex::Real zf = fcx(i,j,k,1);
208#endif
209
210 AMREX_D_TERM(amrex::Real xc = ccc(i,j,k,0);, // centroid of cell (i,j,k)
211 amrex::Real yc = ccc(i,j,k,1);,
212 amrex::Real zc = ccc(i,j,k,2););
213
214 AMREX_D_TERM(amrex::Real delta_x = amrex::Real(0.5) + xc;,
215 amrex::Real delta_y = yf - yc;,
216 amrex::Real delta_z = zf - zc;);
217
218 amrex::Real cc_qmax = amrex::max(q(i,j,k,n),q(i-1,j,k,n));
219 amrex::Real cc_qmin = amrex::min(q(i,j,k,n),q(i-1,j,k,n));
220
221 // Compute slopes of component "n" of q
222 const auto& slopes_eb_hi = amrex_lim_slopes_eb(i, j, k, n, q, ccc, vfrac,
223 AMREX_D_DECL(fcx,fcy,fcz),
224 flag, order);
225
226#if (AMREX_SPACEDIM==3)
227 amrex::Real qpls = q(i ,j,k,n) - delta_x * slopes_eb_hi[0]
228 + delta_y * slopes_eb_hi[1]
229 + delta_z * slopes_eb_hi[2];
230#else
231 amrex::Real qpls = q(i ,j,k,n) - delta_x * slopes_eb_hi[0]
232 + delta_y * slopes_eb_hi[1];
233#endif
234
235 qpls = amrex::max(amrex::min(qpls, cc_qmax), cc_qmin);
236
237 AMREX_D_TERM(xc = ccc(i-1,j,k,0);, // centroid of cell (i-1,j,k)
238 yc = ccc(i-1,j,k,1);,
239 zc = ccc(i-1,j,k,2););
240
241 AMREX_D_TERM(delta_x = amrex::Real(0.5) - xc;,
242 delta_y = yf - yc;,
243 delta_z = zf - zc;);
244
245 // Compute slopes of component "n" of q
246 const auto& slopes_eb_lo = amrex_lim_slopes_eb(i-1, j, k, n, q, ccc, vfrac,
247 AMREX_D_DECL(fcx,fcy,fcz),
248 flag, order);
249
250#if (AMREX_SPACEDIM==3)
251 amrex::Real qmns = q(i-1,j,k,n) + delta_x * slopes_eb_lo[0]
252 + delta_y * slopes_eb_lo[1]
253 + delta_z * slopes_eb_lo[2];
254#else
255 amrex::Real qmns = q(i-1,j,k,n) + delta_x * slopes_eb_lo[0]
256 + delta_y * slopes_eb_lo[1];
257#endif
258
259 qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
260
261 HydroBC::SetEdgeBCsLo(0, i, j, k, n, q, qmns, qpls, umac(i,j,k), d_bcrec[n].lo(0), domain_ilo, is_velocity);
262 HydroBC::SetEdgeBCsHi(0, i, j, k, n, q, qmns, qpls, umac(i,j,k), d_bcrec[n].hi(0), domain_ihi, is_velocity);
263
264 // The clamp below forces the normal velocity to be outflowing at an outflow
265 // face. allow_inflow_on_outflow lets the caller keep an inflowing MAC
266 // velocity there, so the state must not be clamped against it.
267 if (!allow_inflow_on_outflow)
268 {
269 if ( (i==domain_ilo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) )
270 {
271 if ( umac(i,j,k) >= amrex::Real(0.) && n==XVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0.0));
272 qmns = qpls;
273 }
274 if ( (i==domain_ihi+1) && (d_bcrec[n].hi(0) == amrex::BCType::foextrap || d_bcrec[n].hi(0) == amrex::BCType::hoextrap) )
275 {
276 if ( umac(i,j,k) <= amrex::Real(0.) && n==XVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0.0));
277 qpls = qmns;
278 }
279 }
280
281 amrex::Real qs;
282
283 if (umac(i,j,k) > small_vel)
284 {
285 qs = qmns;
286 }
287 else if (umac(i,j,k) < - small_vel)
288 {
289 qs = qpls;
290 }
291 else
292 {
293 qs = amrex::Real(0.5)*(qmns+qpls);
294 }
295
296 return qs;
297}
298
299AMREX_GPU_DEVICE AMREX_FORCE_INLINE
300amrex::Real hydro_ebmol_yedge_state_extdir ( AMREX_D_DECL(int i, int j, int k), int n,
301 amrex::Array4<amrex::Real const> const& q,
302 amrex::Array4<amrex::Real const> const& vmac,
303 AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
304 amrex::Array4<amrex::Real const> const& fcy,
305 amrex::Array4<amrex::Real const> const& fcz),
306 amrex::Array4<amrex::Real const> const& ccc,
307 amrex::Array4<amrex::Real const> const& vfrac,
308 amrex::Array4<amrex::EBCellFlag const> const& flag,
309 amrex::BCRec const* const d_bcrec,
310 amrex::Box const& domain,
311 int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
312{
313
314#if (AMREX_SPACEDIM==2)
315 const int k = 0;
316#endif
317
318 amrex::Real qs;
319 int domlo = domain.smallEnd(1);
320 int domhi = domain.bigEnd(1);
321 bool extlo = d_bcrec[n].lo(1) == amrex::BCType::ext_dir;
322 bool exthi = d_bcrec[n].hi(1) == amrex::BCType::ext_dir;
323
324 // Only short-circuit on faces outside the domain. On the domain face itself we
325 // fall through to SetEdgeBCsLo/Hi, which keeps the Dirichlet value for the
326 // inflowing state but lets the interior state be transported out, so that an
327 // ext_dir face with local backflow upwinds the same way (EB)Godunov does.
328 if ( extlo && j < domlo)
329 {
330 qs = q(i,domlo-1,k,n);
331 }
332 else if (exthi && j > domhi+1)
333 {
334 qs = q(i,domhi+1,k,n);
335 }
336 else
337 {
338 const int domain_ilo = domain.smallEnd(0);
339 const int domain_ihi = domain.bigEnd(0);
340 const int domain_jlo = domain.smallEnd(1);
341 const int domain_jhi = domain.bigEnd(1);
342#if (AMREX_SPACEDIM == 3)
343 const int domain_klo = domain.smallEnd(2);
344 const int domain_khi = domain.bigEnd(2);
345#endif
346
347 AMREX_D_TERM(bool extdir_or_ho_ilo = (d_bcrec[n].lo(0) == amrex::BCType::ext_dir) ||
348 (d_bcrec[n].lo(0) == amrex::BCType::hoextrap);,
349 bool extdir_or_ho_jlo = (d_bcrec[n].lo(1) == amrex::BCType::ext_dir) ||
350 (d_bcrec[n].lo(1) == amrex::BCType::hoextrap);,
351 bool extdir_or_ho_klo = (d_bcrec[n].lo(2) == amrex::BCType::ext_dir) ||
352 (d_bcrec[n].lo(2) == amrex::BCType::hoextrap););
353
354 AMREX_D_TERM(bool extdir_or_ho_ihi = (d_bcrec[n].hi(0) == amrex::BCType::ext_dir) ||
355 (d_bcrec[n].hi(0) == amrex::BCType::hoextrap);,
356 bool extdir_or_ho_jhi = (d_bcrec[n].hi(1) == amrex::BCType::ext_dir) ||
357 (d_bcrec[n].hi(1) == amrex::BCType::hoextrap);,
358 bool extdir_or_ho_khi = (d_bcrec[n].hi(2) == amrex::BCType::ext_dir) ||
359 (d_bcrec[n].hi(2) == amrex::BCType::hoextrap););
360
361 amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of z-face we are extrapolating to
362#if (AMREX_SPACEDIM==3)
363 amrex::Real zf = fcy(i,j,k,1);
364#endif
365
366 AMREX_D_TERM(amrex::Real xc = ccc(i,j,k,0);, // centroid of cell (i,j,k)
367 amrex::Real yc = ccc(i,j,k,1);,
368 amrex::Real zc = ccc(i,j,k,2););
369
370 AMREX_D_TERM(amrex::Real delta_x = xf - xc;,
371 amrex::Real delta_y = amrex::Real(0.5) + yc;,
372 amrex::Real delta_z = zf - zc;);
373
374 amrex::Real cc_qmax = amrex::max(q(i,j,k,n),q(i,j-1,k,n));
375 amrex::Real cc_qmin = amrex::min(q(i,j,k,n),q(i,j-1,k,n));
376
377 // Compute slopes of component "n" of q
378 const auto& slopes_eb_hi =
379 amrex_lim_slopes_extdir_eb(i, j, k, n, q, ccc, vfrac,
380 AMREX_D_DECL(fcx,fcy,fcz), flag,
381 AMREX_D_DECL(extdir_or_ho_ilo, extdir_or_ho_jlo, extdir_or_ho_klo),
382 AMREX_D_DECL(extdir_or_ho_ihi, extdir_or_ho_jhi, extdir_or_ho_khi),
383 AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
384 AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
385 order);
386
387#if (AMREX_SPACEDIM==3)
388 amrex::Real qpls = q(i ,j,k,n) + delta_x * slopes_eb_hi[0]
389 - delta_y * slopes_eb_hi[1]
390 + delta_z * slopes_eb_hi[2];
391#else
392 amrex::Real qpls = q(i ,j,k,n) + delta_x * slopes_eb_hi[0]
393 - delta_y * slopes_eb_hi[1];
394#endif
395
396 qpls = amrex::max(amrex::min(qpls, cc_qmax), cc_qmin);
397
398 AMREX_D_TERM(xc = ccc(i,j-1,k,0);, // centroid of cell (i,j-1,k)
399 yc = ccc(i,j-1,k,1);,
400 zc = ccc(i,j-1,k,2););
401
402 AMREX_D_TERM(delta_x = xf - xc;,
403 delta_y = amrex::Real(0.5) - yc;,
404 delta_z = zf - zc;);
405
406 // Compute slopes of component "n" of q
407 const auto& slopes_eb_lo =
408 amrex_lim_slopes_extdir_eb(i, j-1, k, n, q, ccc, vfrac,
409 AMREX_D_DECL(fcx,fcy,fcz), flag,
410 AMREX_D_DECL(extdir_or_ho_ilo, extdir_or_ho_jlo, extdir_or_ho_klo),
411 AMREX_D_DECL(extdir_or_ho_ihi, extdir_or_ho_jhi, extdir_or_ho_khi),
412 AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
413 AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
414 order);
415
416
417#if (AMREX_SPACEDIM==3)
418 amrex::Real qmns = q(i,j-1,k,n) + delta_x * slopes_eb_lo[0]
419 + delta_y * slopes_eb_lo[1]
420 + delta_z * slopes_eb_lo[2];
421#else
422 amrex::Real qmns = q(i,j-1,k,n) + delta_x * slopes_eb_lo[0]
423 + delta_y * slopes_eb_lo[1];
424#endif
425
426 qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
427
428 HydroBC::SetEdgeBCsLo(1, i, j, k, n, q, qmns, qpls, vmac(i,j,k), d_bcrec[n].lo(1), domain_jlo, is_velocity);
429 HydroBC::SetEdgeBCsHi(1, i, j, k, n, q, qmns, qpls, vmac(i,j,k), d_bcrec[n].hi(1), domain_jhi, is_velocity);
430
431 // The clamp below forces the normal velocity to be outflowing at an outflow
432 // face. allow_inflow_on_outflow lets the caller keep an inflowing MAC
433 // velocity there, so the state must not be clamped against it.
434 if (!allow_inflow_on_outflow)
435 {
436 if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
437 {
438 if ( vmac(i,j,k) >= amrex::Real(0.) && n==YVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0.0));
439 qmns = qpls;
440 }
441 if ( (j==domain_jhi+1) && (d_bcrec[n].hi(1) == amrex::BCType::foextrap || d_bcrec[n].hi(1) == amrex::BCType::hoextrap) )
442 {
443 if ( vmac(i,j,k) <= amrex::Real(0.) && n==YVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0.0));
444 qpls = qmns;
445 }
446 }
447
448 if (vmac(i,j,k) > small_vel)
449 {
450 qs = qmns;
451 }
452 else if (vmac(i,j,k) < - small_vel)
453 {
454 qs = qpls;
455 }
456 else
457 {
458 qs = amrex::Real(0.5)*(qmns+qpls);
459 }
460 }
461
462 return qs;
463}
464
465
466
467AMREX_GPU_DEVICE AMREX_FORCE_INLINE
468amrex::Real hydro_ebmol_yedge_state ( AMREX_D_DECL(int i, int j, int k), int n,
469 amrex::Array4<amrex::Real const> const& q,
470 amrex::Array4<amrex::Real const> const& vmac,
471 AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
472 amrex::Array4<amrex::Real const> const& fcy,
473 amrex::Array4<amrex::Real const> const& fcz),
474 amrex::Array4<amrex::Real const> const& ccc,
475 amrex::Array4<amrex::Real const> const& vfrac,
476 amrex::Array4<amrex::EBCellFlag const> const& flag,
477 amrex::BCRec const* const d_bcrec,
478 amrex::Box const& domain,
479 int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
480{
481#if (AMREX_SPACEDIM==2)
482 const int k = 0;
483#endif
484
485 const int domain_jlo = domain.smallEnd(1);
486 const int domain_jhi = domain.bigEnd(1);
487
488 // local (x,z) of centroid of z-face we are extrapolating to
489 amrex::Real xf = fcy(i,j,k,0);
490#if (AMREX_SPACEDIM==3)
491 amrex::Real zf = fcy(i,j,k,1);
492#endif
493
494 AMREX_D_TERM(amrex::Real xc = ccc(i,j,k,0);, // centroid of cell (i,j,k)
495 amrex::Real yc = ccc(i,j,k,1);,
496 amrex::Real zc = ccc(i,j,k,2););
497
498 AMREX_D_TERM(amrex::Real delta_x = xf - xc;,
499 amrex::Real delta_y = amrex::Real(0.5) + yc;,
500 amrex::Real delta_z = zf - zc;);
501
502 amrex::Real cc_qmax = amrex::max(q(i,j,k,n),q(i,j-1,k,n));
503 amrex::Real cc_qmin = amrex::min(q(i,j,k,n),q(i,j-1,k,n));
504
505 // Compute slopes of component "n" of q
506 const auto& slopes_eb_hi = amrex_lim_slopes_eb(i, j, k, n, q, ccc, vfrac,
507 AMREX_D_DECL(fcx,fcy,fcz),
508 flag, order);
509
510#if (AMREX_SPACEDIM==3)
511 amrex::Real qpls = q(i ,j,k,n) + delta_x * slopes_eb_hi[0]
512 - delta_y * slopes_eb_hi[1]
513 + delta_z * slopes_eb_hi[2];
514#else
515 amrex::Real qpls = q(i ,j,k,n) + delta_x * slopes_eb_hi[0]
516 - delta_y * slopes_eb_hi[1];
517#endif
518
519 qpls = amrex::max(amrex::min(qpls, cc_qmax), cc_qmin);
520
521 AMREX_D_TERM(xc = ccc(i,j-1,k,0);, // centroid of cell (i,j-1,k)
522 yc = ccc(i,j-1,k,1);,
523 zc = ccc(i,j-1,k,2););
524
525 AMREX_D_TERM(delta_x = xf - xc;,
526 delta_y = amrex::Real(0.5) - yc;,
527 delta_z = zf - zc;);
528
529 // Compute slopes of component "n" of q
530 const auto& slopes_eb_lo = amrex_lim_slopes_eb(i, j-1, k, n, q, ccc, vfrac,
531 AMREX_D_DECL(fcx,fcy,fcz),
532 flag, order);
533
534#if (AMREX_SPACEDIM==3)
535 amrex::Real qmns = q(i,j-1,k,n) + delta_x * slopes_eb_lo[0]
536 + delta_y * slopes_eb_lo[1]
537 + delta_z * slopes_eb_lo[2];
538#else
539 amrex::Real qmns = q(i,j-1,k,n) + delta_x * slopes_eb_lo[0]
540 + delta_y * slopes_eb_lo[1];
541#endif
542
543 qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
544
545 HydroBC::SetEdgeBCsLo(1, i, j, k, n, q, qmns, qpls, vmac(i,j,k), d_bcrec[n].lo(1), domain_jlo, is_velocity);
546 HydroBC::SetEdgeBCsHi(1, i, j, k, n, q, qmns, qpls, vmac(i,j,k), d_bcrec[n].hi(1), domain_jhi, is_velocity);
547
548 // The clamp below forces the normal velocity to be outflowing at an outflow
549 // face. allow_inflow_on_outflow lets the caller keep an inflowing MAC
550 // velocity there, so the state must not be clamped against it.
551 if (!allow_inflow_on_outflow)
552 {
553 if ( (j==domain_jlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) )
554 {
555 if ( vmac(i,j,k) >= amrex::Real(0.) && n==YVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0.0));
556 qmns = qpls;
557 }
558 if ( (j==domain_jhi+1) && (d_bcrec[n].hi(1) == amrex::BCType::foextrap || d_bcrec[n].hi(1) == amrex::BCType::hoextrap) )
559 {
560 if ( vmac(i,j,k) <= amrex::Real(0.) && n==YVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0.0));
561 qpls = qmns;
562 }
563 }
564
565 amrex::Real qs;
566
567 if (vmac(i,j,k) > small_vel)
568 {
569 qs = qmns;
570 }
571 else if (vmac(i,j,k) < - small_vel)
572 {
573 qs = qpls;
574 }
575 else
576 {
577 qs = amrex::Real(0.5)*(qmns+qpls);
578 }
579
580 return qs;
581}
582
583#if (AMREX_SPACEDIM==3)
584
585AMREX_GPU_DEVICE AMREX_FORCE_INLINE
586amrex::Real hydro_ebmol_zedge_state_extdir ( int i, int j, int k, int n,
587 amrex::Array4<amrex::Real const> const& q,
588 amrex::Array4<amrex::Real const> const& wmac,
589 AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
590 amrex::Array4<amrex::Real const> const& fcy,
591 amrex::Array4<amrex::Real const> const& fcz),
592 amrex::Array4<amrex::Real const> const& ccc,
593 amrex::Array4<amrex::Real const> const& vfrac,
594 amrex::Array4<amrex::EBCellFlag const> const& flag,
595 amrex::BCRec const* const d_bcrec,
596 amrex::Box const& domain,
597 int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
598{
599 amrex::Real qs;
600 int domlo = domain.smallEnd(2);
601 int domhi = domain.bigEnd(2);
602 bool extlo = d_bcrec[n].lo(2) == amrex::BCType::ext_dir;
603 bool exthi = d_bcrec[n].hi(2) == amrex::BCType::ext_dir;
604
605 // Only short-circuit on faces outside the domain. On the domain face itself we
606 // fall through to SetEdgeBCsLo/Hi, which keeps the Dirichlet value for the
607 // inflowing state but lets the interior state be transported out, so that an
608 // ext_dir face with local backflow upwinds the same way (EB)Godunov does.
609 if ( extlo && k < domlo)
610 {
611 qs = q(i,j,domlo-1,n);
612 }
613 else if ( exthi && k > domhi+1)
614 {
615 qs = q(i,j,domhi+1,n);
616 }
617 else
618 {
619 const int domain_ilo = domain.smallEnd(0);
620 const int domain_ihi = domain.bigEnd(0);
621 const int domain_jlo = domain.smallEnd(1);
622 const int domain_jhi = domain.bigEnd(1);
623#if (AMREX_SPACEDIM == 3)
624 const int domain_klo = domain.smallEnd(2);
625 const int domain_khi = domain.bigEnd(2);
626#endif
627
628 AMREX_D_TERM(bool extdir_or_ho_ilo = (d_bcrec[n].lo(0) == amrex::BCType::ext_dir) ||
629 (d_bcrec[n].lo(0) == amrex::BCType::hoextrap);,
630 bool extdir_or_ho_jlo = (d_bcrec[n].lo(1) == amrex::BCType::ext_dir) ||
631 (d_bcrec[n].lo(1) == amrex::BCType::hoextrap);,
632 bool extdir_or_ho_klo = (d_bcrec[n].lo(2) == amrex::BCType::ext_dir) ||
633 (d_bcrec[n].lo(2) == amrex::BCType::hoextrap););
634
635 AMREX_D_TERM(bool extdir_or_ho_ihi = (d_bcrec[n].hi(0) == amrex::BCType::ext_dir) ||
636 (d_bcrec[n].hi(0) == amrex::BCType::hoextrap);,
637 bool extdir_or_ho_jhi = (d_bcrec[n].hi(1) == amrex::BCType::ext_dir) ||
638 (d_bcrec[n].hi(1) == amrex::BCType::hoextrap);,
639 bool extdir_or_ho_khi = (d_bcrec[n].hi(2) == amrex::BCType::ext_dir) ||
640 (d_bcrec[n].hi(2) == amrex::BCType::hoextrap););
641
642 amrex::Real xf = fcz(i,j,k,0); // local (x,y) of centroid of z-face we are extrapolating to
643 amrex::Real yf = fcz(i,j,k,1);
644
645 amrex::Real xc = ccc(i,j,k,0); // centroid of cell (i,j,k)
646 amrex::Real yc = ccc(i,j,k,1);
647 amrex::Real zc = ccc(i,j,k,2);
648
649 amrex::Real delta_x = xf - xc;
650 amrex::Real delta_y = yf - yc;
651 amrex::Real delta_z = amrex::Real(0.5) + zc;
652
653 amrex::Real cc_qmax = amrex::max(q(i,j,k,n),q(i,j,k-1,n));
654 amrex::Real cc_qmin = amrex::min(q(i,j,k,n),q(i,j,k-1,n));
655
656 // Compute slopes of component "n" of q
657 const auto& slopes_eb_hi =
658 amrex_lim_slopes_extdir_eb(i, j, k, n, q, ccc, vfrac,
659 AMREX_D_DECL(fcx,fcy,fcz), flag,
660 AMREX_D_DECL(extdir_or_ho_ilo, extdir_or_ho_jlo, extdir_or_ho_klo),
661 AMREX_D_DECL(extdir_or_ho_ihi, extdir_or_ho_jhi, extdir_or_ho_khi),
662 AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
663 AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
664 order);
665
666 amrex::Real qpls = q(i,j,k ,n) + delta_x * slopes_eb_hi[0]
667 + delta_y * slopes_eb_hi[1]
668 - delta_z * slopes_eb_hi[2];
669
670 qpls = amrex::max(amrex::min(qpls, cc_qmax), cc_qmin);
671
672 xc = ccc(i,j,k-1,0); // centroid of cell (i,j,k-1)
673 yc = ccc(i,j,k-1,1);
674 zc = ccc(i,j,k-1,2);
675
676 delta_x = xf - xc;
677 delta_y = yf - yc;
678 delta_z = amrex::Real(0.5) - zc;
679
680 // Compute slopes of component "n" of q
681 const auto& slopes_eb_lo =
682 amrex_lim_slopes_extdir_eb(i, j, k-1, n, q, ccc, vfrac,
683 AMREX_D_DECL(fcx,fcy,fcz), flag,
684 AMREX_D_DECL(extdir_or_ho_ilo, extdir_or_ho_jlo, extdir_or_ho_klo),
685 AMREX_D_DECL(extdir_or_ho_ihi, extdir_or_ho_jhi, extdir_or_ho_khi),
686 AMREX_D_DECL(domain_ilo, domain_jlo, domain_klo),
687 AMREX_D_DECL(domain_ihi, domain_jhi, domain_khi),
688 order);
689
690 amrex::Real qmns = q(i,j,k-1,n) + delta_x * slopes_eb_lo[0]
691 + delta_y * slopes_eb_lo[1]
692 + delta_z * slopes_eb_lo[2];
693
694 qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
695
696 HydroBC::SetEdgeBCsLo(2, i, j, k, n, q, qmns, qpls, wmac(i,j,k), d_bcrec[n].lo(2), domain_klo, is_velocity);
697 HydroBC::SetEdgeBCsHi(2, i, j, k, n, q, qmns, qpls, wmac(i,j,k), d_bcrec[n].hi(2), domain_khi, is_velocity);
698
699 // The clamp below forces the normal velocity to be outflowing at an outflow
700 // face. allow_inflow_on_outflow lets the caller keep an inflowing MAC
701 // velocity there, so the state must not be clamped against it.
702 if (!allow_inflow_on_outflow)
703 {
704 if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
705 {
706 if ( wmac(i,j,k) >= amrex::Real(0.) && n==ZVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0.0));
707 qmns = qpls;
708 }
709 if ( (k==domain_khi+1) && (d_bcrec[n].hi(2) == amrex::BCType::foextrap || d_bcrec[n].hi(2) == amrex::BCType::hoextrap) )
710 {
711 if ( wmac(i,j,k) <= amrex::Real(0.) && n==ZVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0.0));
712 qpls = qmns;
713 }
714 }
715
716 if (wmac(i,j,k) > small_vel)
717 {
718 qs = qmns;
719 }
720 else if (wmac(i,j,k) < -small_vel)
721 {
722 qs = qpls;
723 }
724 else
725 {
726 qs = amrex::Real(0.5)*(qmns+qpls);
727 }
728 }
729
730 return qs;
731}
732
733
734
735
736AMREX_GPU_DEVICE AMREX_FORCE_INLINE
737amrex::Real hydro_ebmol_zedge_state ( int i, int j, int k, int n,
738 amrex::Array4<amrex::Real const> const& q,
739 amrex::Array4<amrex::Real const> const& wmac,
740 AMREX_D_DECL(amrex::Array4<amrex::Real const> const& fcx,
741 amrex::Array4<amrex::Real const> const& fcy,
742 amrex::Array4<amrex::Real const> const& fcz),
743 amrex::Array4<amrex::Real const> const& ccc,
744 amrex::Array4<amrex::Real const> const& vfrac,
745 amrex::Array4<amrex::EBCellFlag const> const& flag,
746 amrex::BCRec const* const d_bcrec,
747 amrex::Box const& domain,
748 int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
749{
750 const int domain_klo = domain.smallEnd(2);
751 const int domain_khi = domain.bigEnd(2);
752
753 amrex::Real xf = fcz(i,j,k,0); // local (x,y) of centroid of z-face we are extrapolating to
754 amrex::Real yf = fcz(i,j,k,1);
755
756 amrex::Real xc = ccc(i,j,k,0); // centroid of cell (i,j,k)
757 amrex::Real yc = ccc(i,j,k,1);
758 amrex::Real zc = ccc(i,j,k,2);
759
760 amrex::Real delta_x = xf - xc;
761 amrex::Real delta_y = yf - yc;
762 amrex::Real delta_z = amrex::Real(0.5) + zc;
763
764 amrex::Real cc_qmax = amrex::max(q(i,j,k,n),q(i,j,k-1,n));
765 amrex::Real cc_qmin = amrex::min(q(i,j,k,n),q(i,j,k-1,n));
766
767 // Compute slopes of component "n" of q
768 const auto& slopes_eb_hi = amrex_lim_slopes_eb(i, j, k, n, q, ccc, vfrac,
769 AMREX_D_DECL(fcx,fcy,fcz),
770 flag, order);
771
772 amrex::Real qpls = q(i,j,k ,n) + delta_x * slopes_eb_hi[0]
773 + delta_y * slopes_eb_hi[1]
774 - delta_z * slopes_eb_hi[2];
775
776 qpls = amrex::max(amrex::min(qpls, cc_qmax), cc_qmin);
777
778 xc = ccc(i,j,k-1,0); // centroid of cell (i,j,k-1)
779 yc = ccc(i,j,k-1,1);
780 zc = ccc(i,j,k-1,2);
781
782 delta_x = xf - xc;
783 delta_y = yf - yc;
784 delta_z = amrex::Real(0.5) - zc;
785
786 // Compute slopes of component "n" of q
787 const auto& slopes_eb_lo = amrex_lim_slopes_eb(i, j, k-1, n, q, ccc, vfrac,
788 AMREX_D_DECL(fcx,fcy,fcz),
789 flag, order);
790
791 amrex::Real qmns = q(i,j,k-1,n) + delta_x * slopes_eb_lo[0]
792 + delta_y * slopes_eb_lo[1]
793 + delta_z * slopes_eb_lo[2];
794
795 qmns = amrex::max(amrex::min(qmns, cc_qmax), cc_qmin);
796
797 HydroBC::SetEdgeBCsLo(2, i, j, k, n, q, qmns, qpls, wmac(i,j,k), d_bcrec[n].lo(2), domain_klo, is_velocity);
798 HydroBC::SetEdgeBCsHi(2, i, j, k, n, q, qmns, qpls, wmac(i,j,k), d_bcrec[n].hi(2), domain_khi, is_velocity);
799
800 // The clamp below forces the normal velocity to be outflowing at an outflow
801 // face. allow_inflow_on_outflow lets the caller keep an inflowing MAC
802 // velocity there, so the state must not be clamped against it.
803 if (!allow_inflow_on_outflow)
804 {
805 if ( (k==domain_klo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) )
806 {
807 if ( wmac(i,j,k) >= amrex::Real(0.) && n==ZVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0.0));
808 qmns = qpls;
809 }
810 if ( (k==domain_khi+1) && (d_bcrec[n].hi(2) == amrex::BCType::foextrap || d_bcrec[n].hi(2) == amrex::BCType::hoextrap) )
811 {
812 if ( wmac(i,j,k) <= amrex::Real(0.) && n==ZVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0.0));
813 qpls = qmns;
814 }
815 }
816
817 amrex::Real qs;
818
819 if (wmac(i,j,k) > small_vel)
820 {
821 qs = qmns;
822 }
823 else if (wmac(i,j,k) < -small_vel)
824 {
825 qs = qpls;
826 }
827 else
828 {
829 qs = amrex::Real(0.5)*(qmns+qpls);
830 }
831
832
833 return qs;
834}
835
836#endif
837
838}
839
840#endif
841/** @} */
AMReX_BCRec.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
AMREX_D_TERM
#define AMREX_D_TERM(a, b, c)
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
YVEL
#define YVEL
Definition hydro_constants.H:29
small_vel
static constexpr amrex::Real small_vel
Definition hydro_constants.H:37
XVEL
#define XVEL
Definition hydro_constants.H:28
ZVEL
#define ZVEL
Definition hydro_constants.H:30
amrex::Real
amrex_real Real
amrex::min
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
amrex::max
__host__ __device__ constexpr const T & max(const T &a, const T &b) noexcept
hydro_bcs_K.H
hydro_constants.H
EBMOL::hydro_ebmol_yedge_state
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real hydro_ebmol_yedge_state(AMREX_D_DECL(int i, int j, int k), int n, amrex::Array4< amrex::Real const > const &q, amrex::Array4< amrex::Real const > const &vmac, AMREX_D_DECL(amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz), amrex::Array4< amrex::Real const > const &ccc, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::EBCellFlag const > const &flag, amrex::BCRec const *const d_bcrec, amrex::Box const &domain, int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
Definition hydro_ebmol_edge_state_K.H:468
EBMOL::hydro_ebmol_xedge_state_extdir
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real hydro_ebmol_xedge_state_extdir(AMREX_D_DECL(int i, int j, int k), int n, amrex::Array4< amrex::Real const > const &q, amrex::Array4< amrex::Real const > const &umac, AMREX_D_DECL(amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz), amrex::Array4< amrex::Real const > const &ccc, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::EBCellFlag const > const &flag, amrex::BCRec const *const d_bcrec, amrex::Box const &domain, int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
Definition hydro_ebmol_edge_state_K.H:19
EBMOL::hydro_ebmol_xedge_state
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real hydro_ebmol_xedge_state(AMREX_D_DECL(int i, int j, int k), int n, amrex::Array4< amrex::Real const > const &q, amrex::Array4< amrex::Real const > const &umac, AMREX_D_DECL(amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz), amrex::Array4< amrex::Real const > const &ccc, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::EBCellFlag const > const &flag, amrex::BCRec const *const d_bcrec, amrex::Box const &domain, int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
Definition hydro_ebmol_edge_state_K.H:184
EBMOL::hydro_ebmol_yedge_state_extdir
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real hydro_ebmol_yedge_state_extdir(AMREX_D_DECL(int i, int j, int k), int n, amrex::Array4< amrex::Real const > const &q, amrex::Array4< amrex::Real const > const &vmac, AMREX_D_DECL(amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz), amrex::Array4< amrex::Real const > const &ccc, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::EBCellFlag const > const &flag, amrex::BCRec const *const d_bcrec, amrex::Box const &domain, int order, const bool is_velocity, bool allow_inflow_on_outflow) noexcept
Definition hydro_ebmol_edge_state_K.H:300
HydroBC::SetEdgeBCsHi
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetEdgeBCsHi(int edge_dir, int i, int j, int k, int n, const amrex::Array4< const amrex::Real > &s, amrex::Real &lo, amrex::Real &hi, amrex::Real macvel, int bchi, int domhi, bool is_velocity)
Boundary condition effects.
Definition hydro_bcs_K.H:272
HydroBC::SetEdgeBCsLo
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void SetEdgeBCsLo(int edge_dir, int i, int j, int k, int n, const amrex::Array4< const amrex::Real > &s, amrex::Real &lo, amrex::Real &hi, amrex::Real macvel, int bclo, int domlo, bool is_velocity)
Boundary condition effects.
Definition hydro_bcs_K.H:204