Block-Structured AMR Software Framework
AMReX_EB_Slopes_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_HYDRO_EB_SLOPES_3D_K_H_
2 #define AMREX_HYDRO_EB_SLOPES_3D_K_H_
3 
4 #include <AMReX_Config.H>
5 #include <AMReX_Slopes_K.H>
6 
7 #include <AMReX_EBFArrayBox.H>
8 #include <AMReX_EBCellFlag.H>
9 
10 // amrex_calc_slopes_eb calculates the slope in each coordinate direction using a
11 // least squares linear fit to the 26 nearest neighbors, with the function
12 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
13 // where the data is assume to live, are the same as cell centers. Note that this routine
14 // is called either by amrex_calc_slopes_eb or amrex_calc_slopes_extdir_eb; in the former
15 // A is defined with the cell centroids; in the latter, the A values corresponding to values
16 // defined on faces use the face centroid location.
19 amrex_calc_slopes_eb_given_A (int i, int j, int k, int n,
20  amrex::Real A[27][AMREX_SPACEDIM],
22  amrex::Array4<amrex::EBCellFlag const> const& flag) noexcept
23 {
24  constexpr int dim_a = 27;
25  amrex::Real du[dim_a];
26 
27  int ll=0;
28  for(int kk(-1); kk<=1; kk++) {
29  for(int jj(-1); jj<=1; jj++){
30  for(int ii(-1); ii<=1; ii++){
31  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0))
32  {
33  du[ll] = state(i+ii,j+jj,k+kk,n) - state(i,j,k,n);
34  } else {
35  du[ll] = amrex::Real(0.0);
36  }
37  ll++;
38  }}}
39 
40  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
41  amrex::Real Atb[AMREX_SPACEDIM];
42 
43  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
44  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
45  AtA[jj][ii] = amrex::Real(0.0);
46  }
47  Atb[jj] = amrex::Real(0.0);
48  }
49 
50  for(int lc(0); lc < 27; ++lc)
51  {
52  AtA[0][0] += A[lc][0]* A[lc][0];
53  AtA[0][1] += A[lc][0]* A[lc][1];
54  AtA[0][2] += A[lc][0]* A[lc][2];
55  AtA[1][1] += A[lc][1]* A[lc][1];
56  AtA[1][2] += A[lc][1]* A[lc][2];
57  AtA[2][2] += A[lc][2]* A[lc][2];
58 
59  Atb[0] += A[lc][0]*du[lc];
60  Atb[1] += A[lc][1]*du[lc];
61  Atb[2] += A[lc][2]*du[lc];
62  }
63 
64  // Fill in symmetric
65  AtA[1][0] = AtA[0][1];
66  AtA[2][0] = AtA[0][2];
67  AtA[2][1] = AtA[1][2];
68 
69  amrex::Real detAtA =
70  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
71  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
72  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
73 
74  amrex::Real detAtA_x =
75  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
76  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
77  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
78 
79  // Slope at centroid of (i,j,k)
80  amrex::Real xs = detAtA_x / detAtA;
81 
82  amrex::Real detAtA_y =
83  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
84  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
85  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
86 
87  // Slope at centroid of (i,j,k)
88  amrex::Real ys = detAtA_y / detAtA;
89 
90  amrex::Real detAtA_z =
91  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
92  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
93  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
94 
95  // Slope at centroid of (i,j,k)
96  amrex::Real zs = detAtA_z / detAtA;
97 
98  return {xs,ys,zs};
99 }
100 
101 // amrex_calc_slopes_eb calculates the slope in each coordinate direction using a
102 // least squares linear fit to the 124 nearest neighbors, with the function
103 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
104 // where the data is assume to live, are the same as cell centers. Note that this routine
105 // is called either by amrex_calc_slopes_eb or amrex_calc_slopes_extdir_eb; in the former
106 // A is defined with the cell centroids; in the latter, the A values corresponding to values
107 // defined on faces use the face centroid location.
110 amrex_calc_slopes_eb_given_A_grown (int i, int j, int k, int n, int nx, int ny, int nz,
111  amrex::Real A[125][AMREX_SPACEDIM],
113  amrex::Array4<amrex::EBCellFlag const> const& flag) noexcept
114 {
115  AMREX_ASSERT(nx == 1 || nx == 2);
116  AMREX_ASSERT(ny == 1 || ny == 2);
117  AMREX_ASSERT(nz == 1 || nz == 2);
118 
119  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
120  amrex::Real Atb[AMREX_SPACEDIM];
121 
122  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
123  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
124  AtA[jj][ii] = amrex::Real(0.0);
125  }
126  Atb[jj] = amrex::Real(0.0);
127  }
128 
129  // for(int lc(0); lc < dim_a; ++lc)
130  for (int kk=-nz; kk<=nz; ++kk) {
131  for (int jj=-ny; jj<=ny; ++jj) {
132  for (int ii=-nx; ii<=nx; ++ii) {
133  int lc = (kk+nz)*(2*nx+1)*(2*ny+1) + (jj+ny)*(2*nx+1) + ii+nx;
134  AtA[0][0] += A[lc][0]* A[lc][0];
135  AtA[0][1] += A[lc][0]* A[lc][1];
136  AtA[0][2] += A[lc][0]* A[lc][2];
137  AtA[1][1] += A[lc][1]* A[lc][1];
138  AtA[1][2] += A[lc][1]* A[lc][2];
139  AtA[2][2] += A[lc][2]* A[lc][2];
140 
141  auto du = amrex::Real(0.0);
142  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0) &&
143  ii >= -nx && ii <= nx && jj >= -ny && jj <= ny &&
144  kk >= -nz && kk <= nz) {
145  du = state(i+ii,j+jj,k+kk,n) - state(i,j,k,n);
146  }
147 
148  Atb[0] += A[lc][0]*du;
149  Atb[1] += A[lc][1]*du;
150  Atb[2] += A[lc][2]*du;
151  }
152  }
153  }
154 
155  // Fill in symmetric
156  AtA[1][0] = AtA[0][1];
157  AtA[2][0] = AtA[0][2];
158  AtA[2][1] = AtA[1][2];
159 
160  amrex::Real detAtA =
161  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
162  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
163  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
164 
165  amrex::Real detAtA_x =
166  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
167  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
168  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
169 
170  // Slope at centroid of (i,j,k)
171  amrex::Real xs = detAtA_x / detAtA;
172 
173  amrex::Real detAtA_y =
174  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
175  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
176  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
177 
178  // Slope at centroid of (i,j,k)
179  amrex::Real ys = detAtA_y / detAtA;
180 
181  amrex::Real detAtA_z =
182  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
183  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
184  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
185 
186  // Slope at centroid of (i,j,k)
187  amrex::Real zs = detAtA_z / detAtA;
188 
189  return {xs,ys,zs};
190 }
191 
192 //
193 // amrex_overwrite_with_regular_slopes calculates the slope in each coordinate direction
194 // with a standard non-EB slope calculation (that depends on max_order)
195 //
197 void
198 amrex_overwrite_with_regular_slopes(int i, int j, int k, int n,
199  amrex::Real& xslope, amrex::Real& yslope,
200  amrex::Real& zslope,
203  int max_order) noexcept
204 {
205  //
206  // Here we over-write -- if possible -- with a stencil not using the EB stencil
207  //
208  AMREX_ALWAYS_ASSERT( max_order == 0 || max_order == 2 || max_order == 4);
209 
210  // We have enough cells in the x-direction to do 4th order slopes
211  // centered on (i,j,k) with all values at cell centers
212  if (max_order == 4 &&
213  vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i-2,j,k) == 1. &&
214  vfrac(i+1,j,k) == 1. && vfrac(i+2,j,k) == 1.)
215  {
216  int order = 4;
217  xslope = amrex_calc_xslope(i,j,k,n,order,state);
218  }
219  // We have enough cells in the x-direction to do 2nd order slopes
220  // centered on (i,j,k) with all values at cell centers
221  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.)
222  {
223  int order = 2;
224  xslope = amrex_calc_xslope(i,j,k,n,order,state);
225  } else if (max_order == 0) {
226  xslope = 0.;
227  }
228 
229  // We have enough cells in the y-direction to do 4th order slopes
230  // centered on (i,j,k) with all values at cell centers
231  if (max_order == 4 &&
232  vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j-2,k) == 1. &&
233  vfrac(i,j+1,k) == 1. && vfrac(i,j+2,k) == 1.)
234  {
235  int order = 4;
236  yslope = amrex_calc_yslope(i,j,k,n,order,state);
237  }
238  // We have enough cells in the y-direction to do 2nd order slopes
239  // centered on (i,j,k) with all values at cell centers
240  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.)
241  {
242  int order = 2;
243  yslope = amrex_calc_yslope(i,j,k,n,order,state);
244  } else if (max_order == 0) {
245  yslope = 0.;
246  }
247 
248 #if (AMREX_SPACEDIM == 3)
249  // We have enough cells in the z-direction to do 4th order slopes
250  // centered on (i,j,k) with all values at cell centers
251  if (max_order == 4 &&
252  vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k-2) == 1. &&
253  vfrac(i,j,k+1) == 1. && vfrac(i,j,k+2) == 1.)
254  {
255  int order = 4;
256  zslope = amrex_calc_zslope(i,j,k,n,order,state);
257  }
258  // We have enough cells in the z-direction to do 2nd order slopes
259  // centered on (i,j,k) with all values at cell centers
260  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.)
261  {
262  int order = 2;
263  zslope = amrex_calc_zslope(i,j,k,n,order,state);
264  } else if (max_order == 0) {
265  zslope = 0.;
266  }
267 #endif
268 }
269 
270 // amrex_calc_slopes_eb calculates the slope in each coordinate direction using a
271 // 1) standard limited slope if all three cells in the stencil are regular cells
272 // OR
273 // 2) least squares linear fit to the at-most 26 nearest neighbors, with the function
274 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
275 // where the data is assume to live, are the same as cell centers. Note that calc_slopes_eb
276 // defines the matrix A and passes this A to amrex_calc_slopes_eb_given_A.
277 //
278 // This routine assumes that there are no relevant hoextrap/extdir domain boundary conditions for this cell --
279 // it does not test for them so this should not be called if such boundaries might be present
282 amrex_calc_slopes_eb (int i, int j, int k, int n,
287  int max_order) noexcept
288 {
289  constexpr int dim_a = 27;
290  amrex::Real A[dim_a][AMREX_SPACEDIM];
291 
292  int lc=0;
293  for(int kk(-1); kk<=1; kk++) {
294  for(int jj(-1); jj<=1; jj++) {
295  for(int ii(-1); ii<=1; ii++)
296  {
297  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0))
298  {
299  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) - ccent(i,j,k,0);
300  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) - ccent(i,j,k,1);
301  A[lc][2] = amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) - ccent(i,j,k,2);
302  } else {
303  A[lc][0] = amrex::Real(0.0);
304  A[lc][1] = amrex::Real(0.0);
305  A[lc][2] = amrex::Real(0.0);
306  }
307  lc++;
308  }}}
309 
310  //
311  // These slopes use the EB stencil without testing whether it is actually needed
312  //
313  const auto& slopes = amrex_calc_slopes_eb_given_A (i,j,k,n,A,state,flag);
314  amrex::Real xslope = slopes[0];
315  amrex::Real yslope = slopes[1];
316  amrex::Real zslope = slopes[2];
317 
318  // This will over-write the values of xslope, yslope and/or zslope if appropriate
319  amrex_overwrite_with_regular_slopes(i,j,k,n,xslope,yslope,zslope,state,vfrac,max_order);
320 
321  return {xslope,yslope,zslope};
322 }
323 
324 // amrex_calc_slopes_eb_grown calculates the slope in each coordinate direction using a
325 // 1) standard limited slope if all three cells in the stencil are regular cells
326 // OR
327 // 2) least squares linear fit to the at-most 124 nearest neighbors, with the function
328 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
329 // where the data is assume to live, are the same as cell centers. Note that calc_slopes_eb
330 // defines the matrix A and passes this A to amrex_calc_slopes_eb_given_A.
331 //
332 // This routine assumes that there are no relevant hoextrap/extdir domain boundary conditions for this cell --
333 // it does not test for them so this should not be called if such boundaries might be present
336 amrex_calc_slopes_eb_grown (int i, int j, int k, int n, int nx, int ny, int nz,
341  int max_order) noexcept
342 {
343  AMREX_ASSERT(nx == 1 || nx == 2);
344  AMREX_ASSERT(ny == 1 || ny == 2);
345  AMREX_ASSERT(nz == 1 || nz == 2);
346 
347  constexpr int dim_a = 125;
348  amrex::Real A[dim_a][AMREX_SPACEDIM];
349 
350  // Make sure to zero all the entries in A (since the loop below may not cover all 125)
351  int lc=0;
352  for(int kk(-2); kk<=2; kk++) {
353  for(int jj(-2); jj<=2; jj++) {
354  for(int ii(-2); ii<=2; ii++)
355  {
356  A[lc][0] = amrex::Real(0.0);
357  A[lc][1] = amrex::Real(0.0);
358  A[lc][2] = amrex::Real(0.0);
359  lc++;
360  }}}
361 
362  lc=0;
363  for(int kk(-nz); kk<=nz; kk++) {
364  for(int jj(-ny); jj<=ny; jj++) {
365  for(int ii(-nx); ii<=nx; ii++)
366  {
367  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0))
368  {
369  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) - ccent(i,j,k,0);
370  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) - ccent(i,j,k,1);
371  A[lc][2] = amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) - ccent(i,j,k,2);
372  }
373  lc++;
374  }}}
375  //
376  // These slopes use the EB stencil without testing whether it is actually needed
377  //
378  const auto& slopes = amrex_calc_slopes_eb_given_A_grown (i,j,k,n,nx,ny,nz,A,state,flag);
379  amrex::Real xslope = slopes[0];
380  amrex::Real yslope = slopes[1];
381  amrex::Real zslope = slopes[2];
382 
383  // This will over-write the values of xslope, yslope and/or zslope if appropriate
384  amrex_overwrite_with_regular_slopes(i,j,k,n,xslope,yslope,zslope,state,vfrac,max_order);
385 
386  return {xslope,yslope,zslope};
387 }
388 
389 //
390 // amrex_overwrite_with_regular_slopes_extdir calculates the slope in each coordinate direction
391 // with a standard non-EB slope calculation (that depends on max_order)
392 //
394 void
396  amrex::Real& xslope, amrex::Real& yslope,
397  amrex::Real& zslope,
400  bool edlo_x, bool edlo_y, bool edlo_z,
401  bool edhi_x, bool edhi_y, bool edhi_z,
402  int domlo_x, int domlo_y, int domlo_z,
403  int domhi_x, int domhi_y, int domhi_z,
404  int max_order) noexcept
405 {
406  //
407  // Correct only those cells which are affected by extdir but not by EB:
408  //
409 
410  // We have enough cells in the x-direction to do 4th order slopes
411  // centered on (i,j,k) with all values at cell centers
412  if (max_order == 4 &&
413  vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i-2,j,k) == 1. &&
414  vfrac(i+1,j,k) == 1. && vfrac(i+2,j,k) == 1.)
415  {
416  int order = 4;
417  xslope = amrex_calc_xslope_extdir(i,j,k,n,order,state,edlo_x,edhi_x,domlo_x,domhi_x);
418  }
419  // We have enough cells in the x-direction to do 2nd order slopes
420  // centered on (i,j,k) with all values at cell centers
421  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.)
422  {
423  int order = 2;
424  xslope = amrex_calc_xslope_extdir(i,j,k,n,order,state,edlo_x,edhi_x,domlo_x,domhi_x);
425  } else if (max_order == 0) {
426  xslope = 0.;
427  }
428 
429  // We have enough cells in the y-direction to do 4th order slopes
430  // centered on (i,j,k) with all values at cell centers
431  if (max_order == 4 &&
432  vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j-2,k) == 1. &&
433  vfrac(i,j+1,k) == 1. && vfrac(i,j+2,k) == 1.)
434  {
435  int order = 4;
436  yslope = amrex_calc_yslope_extdir(i,j,k,n,order,state,edlo_y,edhi_y,domlo_y,domhi_y);
437  }
438  // We have enough cells in the y-direction to do 2nd order slopes
439  // centered on (i,j,k) with all values at cell centers
440  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.)
441  {
442  int order = 2;
443  yslope = amrex_calc_yslope_extdir(i,j,k,n,order,state,edlo_y,edhi_y,domlo_y,domhi_y);
444  } else if (max_order == 0) {
445  yslope = 0.;
446  }
447 
448  // We have enough cells in the z-direction to do 4th order slopes
449  // centered on (i,j,k) with all values at cell centers
450  if (max_order == 4 &&
451  vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k-2) == 1. &&
452  vfrac(i,j,k+1) == 1. && vfrac(i,j,k+2) == 1.)
453  {
454  int order = 4;
455  zslope = amrex_calc_zslope_extdir(i,j,k,n,order,state,edlo_z,edhi_z,domlo_z,domhi_z);
456  }
457  // We have enough cells in the z-direction to do 2nd order slopes
458  // centered on (i,j,k) with all values at cell centers
459  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.)
460  {
461  int order = 2;
462  zslope = amrex_calc_zslope_extdir(i,j,k,n,order,state,edlo_z,edhi_z,domlo_z,domhi_z);
463  } else if (max_order == 0) {
464  zslope = 0.;
465  }
466 }
467 
468 // amrex_calc_slopes_extdir_eb calculates the slope in each coordinate direction using a
469 // 1) standard limited slope if all three cells in the stencil are regular cells
470 // (this stencil sees the extdir/hoextrap boundary condition if there is one)
471 // OR
472 // 2) least squares linear fit to the at-most 26 nearest neighbors, with the function
473 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
474 // where the data is assume to live, are the same as cell centers.
475 //
476 // All the slopes are returned in one call.
477 //
478 
481 amrex_calc_slopes_extdir_eb (int i, int j, int k, int n,
489  bool edlo_x, bool edlo_y, bool edlo_z,
490  bool edhi_x, bool edhi_y, bool edhi_z,
491  int domlo_x, int domlo_y, int domlo_z,
492  int domhi_x, int domhi_y, int domhi_z,
493  int max_order) noexcept
494 {
495  constexpr int dim_a = 27;
496 
497  auto xslope = amrex::Real(0.0);
498  auto yslope = amrex::Real(0.0);
499  auto zslope = amrex::Real(0.0);
500 
501  // First get EB-aware slope that doesn't know about extdir
502  bool needs_bdry_stencil = (edlo_x && i <= domlo_x) || (edhi_x && i >= domhi_x) ||
503  (edlo_y && j <= domlo_y) || (edhi_y && j >= domhi_y) ||
504  (edlo_z && k <= domlo_z) || (edhi_z && k >= domhi_z) ;
505 
506  //
507  // This call does not have any knowledge of extdir / hoextrap boundary conditions
508  //
509  if (!needs_bdry_stencil)
510  {
511  // This returns slopes calculated with the regular 1-d approach if all cells in the stencil
512  // are regular. If not, it uses the EB-aware least squares approach to fit a linear profile
513  // using the neighboring un-covered cells.
514  const auto& slopes = amrex_calc_slopes_eb (i,j,k,n,state,ccent,vfrac,flag,max_order);
515  return slopes;
516 
517  } else {
518 
519  amrex::Real A[dim_a][AMREX_SPACEDIM];
520 
521  int lc=0;
522  for(int kk(-1); kk<=1; kk++) {
523  for(int jj(-1); jj<=1; jj++) {
524  for(int ii(-1); ii<=1; ii++)
525  {
526  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0))
527  {
528  bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1);
529  bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1);
530 
531  bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1);
532  bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1);
533 
534  bool klo_test = ( edlo_z && (k == domlo_z) && kk == -1);
535  bool khi_test = ( edhi_z && (k == domhi_z) && kk == 1);
536 
537  // These are the default values if no physical boundary
538  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0);
539  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1);
540  A[lc][2] = amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2);
541  // Do corrections for entire x-face
542  if (ilo_test)
543  {
544  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
545  {
546  A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0);
547  A[lc][2] = amrex::Real(kk) + fcx(i ,j+jj,k+kk,1);
548  }
549  A[lc][0] = -amrex::Real(0.5);
550  } else if (ihi_test) {
551 
552  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
553  {
554  A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0);
555  A[lc][2] = amrex::Real(kk) + fcx(i+ii,j+jj,k+kk,1);
556  }
557  A[lc][0] = amrex::Real(0.5);
558  }
559 
560  // Do corrections for entire y-face
561  if (jlo_test)
562  {
563  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
564  {
565  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0);
566  A[lc][2] = amrex::Real(kk) + fcy(i+ii,j ,k+kk,1);
567  }
568  A[lc][1] = -amrex::Real(0.5);
569 
570  } else if (jhi_test) {
571 
572  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
573  {
574  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0);
575  A[lc][2] = amrex::Real(kk) + fcy(i+ii,j+jj,k+kk,1);
576  }
577  A[lc][1] = amrex::Real(0.5);
578  }
579 
580  // Do corrections for entire z-face
581  if (klo_test)
582  {
583  if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
584  {
585  A[lc][0] = amrex::Real(ii) + fcz(i+ii,j+jj,k ,0);
586  A[lc][1] = amrex::Real(jj) + fcz(i+ii,j+jj,k ,1);
587  }
588  A[lc][2] = -amrex::Real(0.5);
589 
590  } else if (khi_test) {
591  if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
592  {
593  A[lc][0] = amrex::Real(ii) + fcz(i+ii,j+jj,k+kk,0);
594  A[lc][1] = amrex::Real(jj) + fcz(i+ii,j+jj,k+kk,1);
595  }
596  A[lc][2] = amrex::Real(0.5);
597  }
598 
599  A[lc][0] -= ccent(i,j,k,0);
600  A[lc][1] -= ccent(i,j,k,1);
601  A[lc][2] -= ccent(i,j,k,2);
602 
603  } else {
604 
605  A[lc][0] = amrex::Real(0.0);
606  A[lc][1] = amrex::Real(0.0);
607  A[lc][2] = amrex::Real(0.0);
608  }
609  lc++;
610  }}}
611 
612  const auto& slopes = amrex_calc_slopes_eb_given_A (i,j,k,n,A,state,flag);
613  xslope = slopes[0];
614  yslope = slopes[1];
615  zslope = slopes[2];
616 
617  // This will over-write the values of xslope and yslope if appropriate
618  amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,zslope,state,vfrac,
619  edlo_x,edlo_y,edlo_z,edhi_x,edhi_y,edhi_z,
620  domlo_x,domlo_y,domlo_z,domhi_x,domhi_y,domhi_z,
621  max_order);
622 
623  } // end of needs_bndry_stencil
624 
625  // Zero out slopes outside of an extdir (or hoextrap) boundary
626  // TODO: is this the right thing to do at a HOEXTRAP boundary??
627  if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) ||
628  (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) ||
629  (edlo_z && k < domlo_z) || (edhi_z && k > domhi_z) )
630  {
631  xslope = 0.; yslope = 0.; zslope = 0.;
632  }
633 
634  return {xslope,yslope,zslope};
635 }
636 
637 // amrex_calc_slopes_extdir_eb_grown calculates the slope in each coordinate direction using a
638 // 1) standard limited slope if all three cells in the stencil are regular cells
639 // (this stencil sees the extdir/hoextrap boundary condition if there is one)
640 // OR
641 // 2) least squares linear fit to the at-most 124 nearest neighbors, with the function
642 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
643 // where the data is assume to live, are the same as cell centers.
644 //
647 amrex_calc_slopes_extdir_eb_grown (int i, int j, int k, int n,
648  int nx, int ny, int nz,
656  bool edlo_x, bool edlo_y, bool edlo_z,
657  bool edhi_x, bool edhi_y, bool edhi_z,
658  int domlo_x, int domlo_y, int domlo_z,
659  int domhi_x, int domhi_y, int domhi_z,
660  int max_order) noexcept
661 {
662  constexpr int dim_a = 125;
663 
664  auto xslope = amrex::Real(0.0);
665  auto yslope = amrex::Real(0.0);
666  auto zslope = amrex::Real(0.0);
667 
668  // First get EB-aware slope that doesn't know about extdir
669  bool needs_bdry_stencil = (edlo_x && i <= domlo_x) || (edhi_x && i >= domhi_x) ||
670  (edlo_y && j <= domlo_y) || (edhi_y && j >= domhi_y) ||
671  (edlo_z && k <= domlo_z) || (edhi_z && k >= domhi_z) ;
672 
673  //
674  // This call does not have any knowledge of extdir / hoextrap boundary conditions
675  //
676  if (!needs_bdry_stencil)
677  {
678  // This returns slopes calculated with the regular 1-d approach if all cells in the stencil
679  // are regular. If not, it uses the EB-aware least squares approach to fit a linear profile
680  // using the neighboring un-covered cells.
681  const auto& slopes = amrex_calc_slopes_eb_grown (i,j,k,n,nx,ny,nz,state,ccent,vfrac,flag,max_order);
682  return slopes;
683 
684  } else {
685 
686  amrex::Real A[dim_a][AMREX_SPACEDIM];
687 
688  int lc=0;
689  for(int kk(-nz); kk<=nz; kk++) {
690  for(int jj(-ny); jj<=ny; jj++) {
691  for(int ii(-nx); ii<=nx; ii++)
692  {
693  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0))
694  {
695  bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1);
696  bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1);
697 
698  bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1);
699  bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1);
700 
701  bool klo_test = ( edlo_z && (k == domlo_z) && kk == -1);
702  bool khi_test = ( edhi_z && (k == domhi_z) && kk == 1);
703 
704  // These are the default values if no physical boundary
705  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0);
706  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1);
707  A[lc][2] = amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2);
708  // Do corrections for entire x-face
709  if (ilo_test)
710  {
711  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
712  {
713  A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0);
714  A[lc][2] = amrex::Real(kk) + fcx(i ,j+jj,k+kk,1);
715  }
716  A[lc][0] = -amrex::Real(0.5);
717  } else if (ihi_test) {
718 
719  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
720  {
721  A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0);
722  A[lc][2] = amrex::Real(kk) + fcx(i+ii,j+jj,k+kk,1);
723  }
724  A[lc][0] = amrex::Real(0.5);
725  }
726 
727  // Do corrections for entire y-face
728  if (jlo_test)
729  {
730  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
731  {
732  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0);
733  A[lc][2] = amrex::Real(kk) + fcy(i+ii,j ,k+kk,1);
734  }
735  A[lc][1] = -amrex::Real(0.5);
736 
737  } else if (jhi_test) {
738 
739  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
740  {
741  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0);
742  A[lc][2] = amrex::Real(kk) + fcy(i+ii,j+jj,k+kk,1);
743  }
744  A[lc][1] = amrex::Real(0.5);
745  }
746 
747  // Do corrections for entire z-face
748  if (klo_test)
749  {
750  if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
751  {
752  A[lc][0] = amrex::Real(ii) + fcz(i+ii,j+jj,k ,0);
753  A[lc][1] = amrex::Real(jj) + fcz(i+ii,j+jj,k ,1);
754  }
755  A[lc][2] = -amrex::Real(0.5);
756 
757  } else if (khi_test) {
758  if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
759  {
760  A[lc][0] = amrex::Real(ii) + fcz(i+ii,j+jj,k+kk,0);
761  A[lc][1] = amrex::Real(jj) + fcz(i+ii,j+jj,k+kk,1);
762  }
763  A[lc][2] = amrex::Real(0.5);
764  }
765 
766  A[lc][0] -= ccent(i,j,k,0);
767  A[lc][1] -= ccent(i,j,k,1);
768  A[lc][2] -= ccent(i,j,k,2);
769 
770  } else {
771  A[lc][0] = amrex::Real(0.0);
772  A[lc][1] = amrex::Real(0.0);
773  A[lc][2] = amrex::Real(0.0);
774  }
775  lc++;
776  }}}
777 
778  const auto& slopes = amrex_calc_slopes_eb_given_A_grown (i,j,k,n,nx,ny,nz,A,state,flag);
779  xslope = slopes[0];
780  yslope = slopes[1];
781  zslope = slopes[2];
782 
783  // This will over-write the values of xslope and yslope if appropriate
784  amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,zslope,state,vfrac,
785  edlo_x,edlo_y,edlo_z,edhi_x,edhi_y,edhi_z,
786  domlo_x,domlo_y,domlo_z,domhi_x,domhi_y,domhi_z,
787  max_order);
788 
789  } // end of needs_bndry_stencil
790 
791  // Zero out slopes outside of an extdir (or hoextrap) boundary
792  // TODO: is this the right thing to do at a HOEXTRAP boundary??
793  if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) ||
794  (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) ||
795  (edlo_z && k < domlo_z) || (edhi_z && k > domhi_z) )
796  {
797  xslope = 0.; yslope = 0.; zslope = 0.;
798  }
799 
800  return {xslope,yslope,zslope};
801 }
802 
804 amrex::Real
805 amrex_calc_alpha_stencil(amrex::Real q_hat, amrex::Real q_max,
806  amrex::Real q_min, amrex::Real state, amrex::Real alpha) noexcept
807 {
808  using namespace amrex::literals;
809 
810  auto alpha_temp = amrex::Real(0.0);
811  auto small = amrex::Real(1.0e-13);
812 
813  if ((q_hat-state) > small) {
814  alpha_temp = amrex::min(1.0_rt,(q_max-state)/(q_hat-state));
815  } else if ((q_hat-state) < -small) {
816  alpha_temp = amrex::min(1.0_rt,(q_min-state)/(q_hat-state));
817  } else {
818  alpha_temp = 1.0_rt;
819  }
820  return amrex::min(alpha, alpha_temp);
821 }
822 
825 amrex_calc_alpha_limiter(int i, int j, int k, int n,
832  amrex::Array4<amrex::Real const> const& ccent) noexcept
833 {
834  auto alpha = amrex::Real(2.0);
835 
836  int cuts_x = 0; int cuts_y = 0; int cuts_z = 0;
837 
838  // Compute how many cut or regular faces do we have in 3x3 block
839  for(int kk(-1); kk<=1; kk++)
840  {
841  for(int jj(-1); jj<=1; jj++){
842  for(int ii(-1); ii<=1; ii++){
843  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0))
844  {
845  if ((ii==-1 || ii==1) && jj==0 && kk==0) { cuts_x++; }
846  if ((jj==-1 || jj==1) && ii==0 && kk==0) { cuts_y++; }
847  if ((kk==-1 || kk==1) && ii==0 && jj==0) { cuts_z++; }
848  }
849  }
850  }
851  }
852 
853  amrex::Real xc = ccent(i,j,k,0); // centroid of cell (i,j,k)
854  amrex::Real yc = ccent(i,j,k,1);
855  amrex::Real zc = ccent(i,j,k,2);
856 
857 
858  //Reconstruct values at the face centroids and compute the limiter
859  if (flag(i,j,k).isConnected(0,1,0)) {
860  amrex::Real xf = fcy(i,j+1,k,0); // local (x,z) of centroid of y-face we are extrapolating to
861  amrex::Real zf = fcy(i,j+1,k,1); // local (x,z) of centroid of y-face we are extrapolating to
862 
863  amrex::Real delta_x = xf - xc;
864  amrex::Real delta_y = amrex::Real(0.5) - yc;
865  amrex::Real delta_z = zf - zc;
866 
867  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
868  + delta_y * slopes[1]
869  + delta_z * slopes[2];
870 
871  // Compute max and min values in a 3x2x3 stencil
872  amrex::Real q_min = state(i,j,k,n);
873  amrex::Real q_max = state(i,j,k,n);
874  for(int kk(-1); kk<=1; kk++) {
875  for(int jj(0); jj<=1; jj++){
876  for(int ii(-1); ii<=1; ii++){
877  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
878  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
879  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
880  }
881  }
882  }
883  }
884 
885  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
886  }
887  if (flag(i,j,k).isConnected(0,-1,0)){
888  amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of y-face we are extrapolating to
889  amrex::Real zf = fcy(i,j,k,1); // local (x,z) of centroid of y-face we are extrapolating to
890 
891  amrex::Real delta_x = xf - xc;
892  amrex::Real delta_y = amrex::Real(0.5) + yc;
893  amrex::Real delta_z = zf - zc;
894 
895  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
896  - delta_y * slopes[1]
897  + delta_z * slopes[2];
898 
899  // Compute max and min values in a 3x2x3 stencil
900  amrex::Real q_min = state(i,j,k,n);
901  amrex::Real q_max = state(i,j,k,n);
902  for(int kk(-1); kk<=1; kk++) {
903  for(int jj(-1); jj<=0; jj++){
904  for(int ii(-1); ii<=1; ii++){
905  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
906  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
907  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
908  }
909  }
910  }
911  }
912 
913  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
914  }
915  if (flag(i,j,k).isConnected(1,0,0)) {
916  amrex::Real yf = fcx(i+1,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to
917  amrex::Real zf = fcx(i+1,j,k,1); // local (y,z) of centroid of x-face we are extrapolating to
918 
919  amrex::Real delta_x = amrex::Real(0.5) - xc;
920  amrex::Real delta_y = yf - yc;
921  amrex::Real delta_z = zf - zc;
922 
923  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
924  + delta_y * slopes[1]
925  + delta_z * slopes[2];
926 
927  // Compute max and min values in a 2x3x3 stencil
928  amrex::Real q_min = state(i,j,k,n);
929  amrex::Real q_max = state(i,j,k,n);
930  for(int kk(-1); kk<=1; kk++) {
931  for(int jj(-1); jj<=1; jj++){
932  for(int ii(0); ii<=1; ii++){
933  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
934  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
935  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
936  }
937  }
938  }
939  }
940 
941  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
942  }
943  if (flag(i,j,k).isConnected(-1,0,0)) {
944  amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to
945  amrex::Real zf = fcx(i,j,k,1); // local (y,z) of centroid of x-face we are extrapolating to
946 
947  amrex::Real delta_x = amrex::Real(0.5) + xc;
948  amrex::Real delta_y = yf - yc;
949  amrex::Real delta_z = zf - zc;
950 
951  amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0]
952  + delta_y * slopes[1]
953  + delta_z * slopes[2];
954 
955  // Compute max and min values in a 3x2x3 stencil
956  amrex::Real q_min = state(i,j,k,n);
957  amrex::Real q_max = state(i,j,k,n);
958  for(int kk(-1); kk<=1; kk++) {
959  for(int jj(-1); jj<=1; jj++){
960  for(int ii(-1); ii<=0; ii++){
961  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
962  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
963  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
964  }
965  }
966  }
967  }
968  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
969  }
970  if (flag(i,j,k).isConnected(0,0,1)) {
971  amrex::Real xf = fcz(i,j,k+1,0); // local (x,y) of centroid of z-face we are extrapolating to
972  amrex::Real yf = fcz(i,j,k+1,1); // local (x,y) of centroid of z-face we are extrapolating to
973 
974  amrex::Real delta_x = xf - xc;
975  amrex::Real delta_y = yf - yc;
976  amrex::Real delta_z = amrex::Real(0.5) - zc;
977 
978  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
979  + delta_y * slopes[1]
980  + delta_z * slopes[2];
981 
982  // Compute max and min values in a 3x3x2 stencil
983  amrex::Real q_min = state(i,j,k,n);
984  amrex::Real q_max = state(i,j,k,n);
985  for(int kk(0); kk<=1; kk++) {
986  for(int jj(-1); jj<=1; jj++){
987  for(int ii(-1); ii<=1; ii++){
988  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
989  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
990  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
991  }
992  }
993  }
994  }
995  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
996  }
997  if (flag(i,j,k).isConnected(0,0,-1)) {
998  amrex::Real xf = fcz(i,j,k,0); // local (x,y) of centroid of z-face we are extrapolating to
999  amrex::Real yf = fcz(i,j,k,1); // local (x,y) of centroid of z-face we are extrapolating to
1000 
1001  amrex::Real delta_x = xf - xc;
1002  amrex::Real delta_y = yf - yc;
1003  amrex::Real delta_z = amrex::Real(0.5) + zc;
1004 
1005  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
1006  + delta_y * slopes[1]
1007  - delta_z * slopes[2];
1008 
1009  // Compute max and min values in a 3x2x3 stencil
1010  amrex::Real q_min = state(i,j,k,n);
1011  amrex::Real q_max = state(i,j,k,n);
1012  for(int kk(-1); kk<=0; kk++) {
1013  for(int jj(-1); jj<=1; jj++){
1014  for(int ii(-1); ii<=1; ii++){
1015  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
1016  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
1017  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
1018  }
1019  }
1020  }
1021  }
1022  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
1023  }
1024 
1025  amrex::Real xalpha = alpha;
1026  amrex::Real yalpha = alpha;
1027  amrex::Real zalpha = alpha;
1028 
1029  //Zeroing out the slopes in the direction where a covered face exists.
1030  if (cuts_x<2) { xalpha = 0; }
1031  if (cuts_y<2) { yalpha = 0; }
1032  if (cuts_z<2) { zalpha = 0; }
1033 
1034  return {xalpha,yalpha,zalpha};
1035 }
1036 
1037 //amrex_lim_slopes_eb computes the slopes calling amrex_calc_slopes_eb, and then each slope component
1038 //is multiplied by a limiter based on the work of Barth-Jespersen.
1041 amrex_lim_slopes_eb (int i, int j, int k, int n,
1042  amrex::Array4<amrex::Real const> const& state,
1043  amrex::Array4<amrex::Real const> const& ccent,
1044  amrex::Array4<amrex::Real const> const& vfrac,
1049  int max_order) noexcept
1050 {
1051 
1054 
1055  slopes = amrex_calc_slopes_eb(i,j,k,n,state,ccent,vfrac,flag,max_order);
1056 
1057  alpha_lim = amrex_calc_alpha_limiter(i,j,k,n,state,flag,slopes,fcx,fcy,fcz,ccent);
1058 
1059  // Setting limiter to 1 for stencils that just consists of non-EB cells because
1060  // amrex_calc_slopes_eb routine will call the slope routine for non-EB stencils that has already a limiter
1061  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) {
1062  alpha_lim[0] = 1.0;
1063  }
1064 
1065  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) {
1066  alpha_lim[1] = 1.0;
1067  }
1068 
1069  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) {
1070  alpha_lim[2] = 1.0;
1071  }
1072 
1073  return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1],alpha_lim[2]*slopes[2]};
1074 }
1075 
1076 //amrex_lim_slopes_extdir_eb computes the slopes calling amrex_calc_slopes_extdir_eb, and then each slope component
1077 //is multiplied by a limiter based on the work of Barth-Jespersen.
1080 amrex_lim_slopes_extdir_eb (int i, int j, int k, int n,
1081  amrex::Array4<amrex::Real const> const& state,
1082  amrex::Array4<amrex::Real const> const& ccent,
1083  amrex::Array4<amrex::Real const> const& vfrac,
1088  bool edlo_x, bool edlo_y, bool edlo_z,
1089  bool edhi_x, bool edhi_y, bool edhi_z,
1090  int domlo_x, int domlo_y, int domlo_z,
1091  int domhi_x, int domhi_y, int domhi_z,
1092  int max_order) noexcept
1093 {
1094 
1097 
1098  slopes = amrex_calc_slopes_extdir_eb(i,j,k,n,state,ccent,vfrac,fcx,fcy,fcz,flag,
1099  edlo_x,edlo_y,edlo_z,edhi_x,edhi_y,edhi_z,
1100  domlo_x,domlo_y,domlo_z,domhi_x,domhi_y,domhi_z,max_order);
1101  alpha_lim = amrex_calc_alpha_limiter(i,j,k,n,state,flag,slopes,fcx,fcy,fcz,ccent);
1102 
1103  // Setting limiter to 1 for stencils that just consists of non-EB cells because
1104  // amrex_calc_slopes_extdir_eb routine will call the slope routine for non-EB stencils that has already a limiter
1105  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) {
1106  alpha_lim[0] = 1.0;
1107  }
1108 
1109  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) {
1110  alpha_lim[1] = 1.0;
1111  }
1112 
1113  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) {
1114  alpha_lim[2] = 1.0;
1115  }
1116 
1117  return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1],alpha_lim[2]*slopes[2]};
1118 }
1119 
1120 #endif
#define AMREX_ASSERT(EX)
Definition: AMReX_BLassert.H:38
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_eb_given_A_grown(int i, int j, int k, int n, int nx, int ny, int nz, amrex::Real A[125][AMREX_SPACEDIM], amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::EBCellFlag const > const &flag) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:110
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void amrex_overwrite_with_regular_slopes_extdir(int i, int j, int k, int n, amrex::Real &xslope, amrex::Real &yslope, amrex::Real &zslope, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &vfrac, bool edlo_x, bool edlo_y, bool edlo_z, bool edhi_x, bool edhi_y, bool edhi_z, int domlo_x, int domlo_y, int domlo_z, int domhi_x, int domhi_y, int domhi_z, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:395
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_lim_slopes_eb(int i, int j, int k, int n, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz, amrex::Array4< amrex::EBCellFlag const > const &flag, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:1041
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_extdir_eb(int i, int j, int k, int n, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz, amrex::Array4< amrex::EBCellFlag const > const &flag, bool edlo_x, bool edlo_y, bool edlo_z, bool edhi_x, bool edhi_y, bool edhi_z, int domlo_x, int domlo_y, int domlo_z, int domhi_x, int domhi_y, int domhi_z, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:481
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real amrex_calc_alpha_stencil(amrex::Real q_hat, amrex::Real q_max, amrex::Real q_min, amrex::Real state, amrex::Real alpha) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:805
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void amrex_overwrite_with_regular_slopes(int i, int j, int k, int n, amrex::Real &xslope, amrex::Real &yslope, amrex::Real &zslope, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &vfrac, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:198
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_eb_grown(int i, int j, int k, int n, int nx, int ny, int nz, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::EBCellFlag const > const &flag, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:336
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_lim_slopes_extdir_eb(int i, int j, int k, int n, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz, amrex::Array4< amrex::EBCellFlag const > const &flag, bool edlo_x, bool edlo_y, bool edlo_z, bool edhi_x, bool edhi_y, bool edhi_z, int domlo_x, int domlo_y, int domlo_z, int domhi_x, int domhi_y, int domhi_z, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:1080
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_extdir_eb_grown(int i, int j, int k, int n, int nx, int ny, int nz, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::Real const > const &fcx, amrex::Array4< amrex::Real const > const &fcy, amrex::Array4< amrex::Real const > const &fcz, amrex::Array4< amrex::EBCellFlag const > const &flag, bool edlo_x, bool edlo_y, bool edlo_z, bool edhi_x, bool edhi_y, bool edhi_z, int domlo_x, int domlo_y, int domlo_z, int domhi_x, int domhi_y, int domhi_z, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:647
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_eb(int i, int j, int k, int n, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::Real const > const &vfrac, amrex::Array4< amrex::EBCellFlag const > const &flag, int max_order) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:282
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_eb_given_A(int i, int j, int k, int n, amrex::Real A[27][AMREX_SPACEDIM], amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::EBCellFlag const > const &flag) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:19
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_alpha_limiter(int i, int j, int k, int n, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::EBCellFlag const > const &flag, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &slopes, 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 &ccent) noexcept
Definition: AMReX_EB_Slopes_3D_K.H:825
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
AMREX_GPU_DEVICE AMREX_FORCE_INLINE Real amrex_calc_yslope_extdir(int i, int j, int k, int n, int order, amrex::Array4< Real const > const &q, bool edlo, bool edhi, int domlo, int domhi) noexcept
Definition: AMReX_Slopes_K.H:200
AMREX_GPU_DEVICE AMREX_FORCE_INLINE Real amrex_calc_yslope(int i, int j, int k, int n, int order, amrex::Array4< Real const > const &q) noexcept
Definition: AMReX_Slopes_K.H:151
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
AMREX_GPU_DEVICE AMREX_FORCE_INLINE Real amrex_calc_xslope_extdir(int i, int j, int k, int n, int order, amrex::Array4< Real const > const &q, bool edlo, bool edhi, int domlo, int domhi) noexcept
Definition: AMReX_Slopes_K.H:60
AMREX_GPU_DEVICE AMREX_FORCE_INLINE Real amrex_calc_xslope(int i, int j, int k, int n, int order, amrex::Array4< Real const > const &q) noexcept
Definition: AMReX_Slopes_K.H:10
Definition: AMReX_Array4.H:61
Definition: AMReX_Array.H:34