Block-Structured AMR Software Framework
AMReX_EB_Slopes_2D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_HYDRO_EB_SLOPES_2D_K_H_
2 #define AMREX_HYDRO_EB_SLOPES_2D_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 8 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.
17 //
20 amrex_calc_slopes_eb_given_A (int i, int j, int /*k*/, int n,
21  amrex::Real A[9][AMREX_SPACEDIM],
23  amrex::Array4<amrex::EBCellFlag const> const& flag) noexcept
24 {
25  constexpr int dim_a = 9;
26  amrex::Real du[dim_a];
27 
28  int ll=0;
29  for(int jj(-1); jj<=1; jj++){
30  for(int ii(-1); ii<=1; ii++){
31  if( flag(i,j,0).isConnected(ii,jj,0) &&
32  ! (ii==0 && jj==0)) {
33  du[ll] = state(i+ii,j+jj,0,n) - state(i,j,0,n);
34  } else {
35  du[ll] = amrex::Real(0.0);
36  }
37  ll++;
38  }
39  }
40 
41  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
42  amrex::Real Atb[AMREX_SPACEDIM];
43 
44  for(int jj(0); jj<AMREX_SPACEDIM; ++jj)
45  {
46  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
47  AtA[jj][ii] = amrex::Real(0.0);
48  }
49  Atb[jj] = amrex::Real(0.0);
50  }
51 
52  for(int lc(0); lc<dim_a; ++lc)
53  {
54  AtA[0][0] += A[lc][0]* A[lc][0];
55  AtA[0][1] += A[lc][0]* A[lc][1];
56  AtA[1][1] += A[lc][1]* A[lc][1];
57 
58  Atb[0] += A[lc][0]*du[lc];
59  Atb[1] += A[lc][1]*du[lc];
60  }
61 
62  // Fill in symmetric
63  AtA[1][0] = AtA[0][1];
64 
65  amrex::Real detAtA =
66  (AtA[0][0]*AtA[1][1])-
67  (AtA[0][1]*AtA[1][0]);
68 
69  amrex::Real detAtA_x =
70  (Atb[0] *AtA[1][1]) -
71  (AtA[0][1]*Atb[1]);
72 
73  // Slope at centroid of (i,j,k)
74  amrex::Real xs = detAtA_x / detAtA;
75 
76  amrex::Real detAtA_y =
77  (AtA[0][0]*Atb[1]) -
78  (Atb[0] * AtA[1][0]);
79 
80  // Slope at centroid of (i,j,k)
81  amrex::Real ys = detAtA_y / detAtA;
82 
83  return {xs,ys};
84 }
85 
86 // amrex_calc_slopes_eb_grown calculates the slope in each coordinate direction using a
87 // least squares linear fit to the (up to) 24 nearest neighbors, with the function
88 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
89 // where the data is assume to live, are the same as cell centers. Note that this routine
90 // is called either by amrex_calc_slopes_eb_grown or amrex_calc_slopes_extdir_eb_grown; in the former
91 // A is defined with the cell centroids; in the latter, the A values corresponding to values
92 // defined on faces use the face centroid location.
93 //
96 amrex_calc_slopes_eb_given_A_grown (int i, int j, int /*k*/, int n, int nx, int ny,
97  amrex::Real A[25][AMREX_SPACEDIM],
99  amrex::Array4<amrex::EBCellFlag const> const& flag) noexcept
100 {
101  AMREX_ASSERT(nx == 1 || nx == 2);
102  AMREX_ASSERT(ny == 1 || ny == 2);
103 
104  constexpr int dim_a = 25;
105 
106  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
107  amrex::Real Atb[AMREX_SPACEDIM];
108 
109  for(int jj(0); jj<AMREX_SPACEDIM; ++jj)
110  {
111  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
112  AtA[jj][ii] = amrex::Real(0.0);
113  }
114  Atb[jj] = amrex::Real(0.0);
115  }
116 
117  amrex::Real du;
118  for(int jj(-ny); jj<=ny; jj++){
119  for(int ii(-nx); ii<=nx; ii++){
120  int lc = (jj+ny)*(2*nx+1) + ii+nx;
121  if(!flag(i+ii,j+jj,0).isCovered() && !(ii==0 && jj==0)) {
122  du = state(i+ii,j+jj,0,n) - state(i,j,0,n);
123  Atb[0] += A[lc][0]*du;
124  Atb[1] += A[lc][1]*du;
125  }
126  lc++;
127  }
128  }
129 
130  for(int lc(0); lc<dim_a; ++lc)
131  {
132  AtA[0][0] += A[lc][0]* A[lc][0];
133  AtA[0][1] += A[lc][0]* A[lc][1];
134  AtA[1][1] += A[lc][1]* A[lc][1];
135  }
136 
137  // Fill in symmetric
138  AtA[1][0] = AtA[0][1];
139 
140  amrex::Real detAtA =
141  (AtA[0][0]*AtA[1][1])-
142  (AtA[0][1]*AtA[1][0]);
143 
144  amrex::Real detAtA_x =
145  (Atb[0] *AtA[1][1]) -
146  (AtA[0][1]*Atb[1]);
147 
148  // Slope at centroid of (i,j,k)
149  amrex::Real xs = detAtA_x / detAtA;
150 
151  amrex::Real detAtA_y =
152  (AtA[0][0]*Atb[1]) -
153  (Atb[0] * AtA[1][0]);
154 
155  // Slope at centroid of (i,j,k)
156  amrex::Real ys = detAtA_y / detAtA;
157 
158  return {xs,ys};
159 }
160 
161 //
162 // amrex_overwrite_with_regular_slopes calculates the slope in each coordinate direction
163 // with a standard non-EB slope calculation (that depends on max_order)
164 //
166 void
167 amrex_overwrite_with_regular_slopes(int i, int j, int k, int n,
168  amrex::Real& xslope, amrex::Real& yslope,
171  int max_order) noexcept
172 {
173  //
174  // Here we over-write -- if possible -- with a stencil not using the EB stencil
175  //
176  AMREX_ALWAYS_ASSERT( max_order == 0 || max_order == 2 || max_order == 4);
177 
178  // We have enough cells in the x-direction to do 4th order slopes
179  // centered on (i,j,k) with all values at cell centers
180  if (max_order == 4 &&
181  vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i-2,j,k) == 1. &&
182  vfrac(i+1,j,k) == 1. && vfrac(i+2,j,k) == 1.)
183  {
184  int order = 4;
185  xslope = amrex_calc_xslope(i,j,k,n,order,state);
186  }
187  // We have enough cells in the x-direction to do 2nd order slopes
188  // centered on (i,j,k) with all values at cell centers
189  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.)
190  {
191  int order = 2;
192  xslope = amrex_calc_xslope(i,j,k,n,order,state);
193  } else if (max_order == 0) {
194  xslope = 0.;
195  }
196 
197  // We have enough cells in the y-direction to do 4th order slopes
198  // centered on (i,j,k) with all values at cell centers
199  if (max_order == 4 &&
200  vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j-2,k) == 1. &&
201  vfrac(i,j+1,k) == 1. && vfrac(i,j+2,k) == 1.)
202  {
203  int order = 4;
204  yslope = amrex_calc_yslope(i,j,k,n,order,state);
205  }
206  // We have enough cells in the y-direction to do 2nd order slopes
207  // centered on (i,j,k) with all values at cell centers
208  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.)
209  {
210  int order = 2;
211  yslope = amrex_calc_yslope(i,j,k,n,order,state);
212  } else if (max_order == 0) {
213  yslope = 0.;
214  }
215 }
216 
217 // amrex_calc_slopes_eb calculates the slope in each coordinate direction using a
218 // 1) standard limited slope if all three cells in the stencil are regular cells
219 // OR
220 // 2) least squares linear fit to the at-most 8 nearest neighbors, with the function
221 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
222 // where the data is assume to live, are the same as cell centers. Note that calc_slopes_eb
223 // defines the matrix A and passes this A to amrex_calc_slopes_eb_given_A.
224 //
225 // This routine assumes that there are no relevant hoextrap/extdir domain boundary conditions for this cell --
226 // it does not test for them so this should not be called if such boundaries might be present
227 //
230 amrex_calc_slopes_eb (int i, int j, int k, int n,
235  int max_order) noexcept
236 {
237  constexpr int dim_a = 9;
238  amrex::Real A[dim_a][AMREX_SPACEDIM];
239 
240  int lc=0;
241  int kk = 0;
242  {
243  for(int jj(-1); jj<=1; jj++){
244  for(int ii(-1); ii<=1; ii++){
245 
246  if( flag(i,j,k).isConnected(ii,jj,kk) &&
247  ! (ii==0 && jj==0 && kk==0)) {
248 
249  // Not multiplying by dx to be consistent with how the
250  // slope is stored. Also not including the global shift
251  // wrt plo or i,j,k. We only need relative distance.
252 
253  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) - ccent(i,j,k,0);
254  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) - ccent(i,j,k,1);
255 
256  } else {
257 
258  A[lc][0] = amrex::Real(0.0);
259  A[lc][1] = amrex::Real(0.0);
260  }
261  lc++;
262  }
263  }
264  }
265 
266  const auto& eb_slopes = amrex_calc_slopes_eb_given_A (i,j,k,n,A,state,flag);
267 
268  amrex::Real xslope = eb_slopes[0];
269  amrex::Real yslope = eb_slopes[1];
270 
271  // This will over-write the values of xslope and yslope if appropriate
272  amrex_overwrite_with_regular_slopes(i,j,k,n,xslope,yslope,state,vfrac,max_order);
273 
274  return {xslope,yslope};
275 }
276 
277 // amrex_calc_slopes_eb_grown calculates the slope in each coordinate direction using a
278 // 1) standard limited slope if all three cells in the stencil are regular cells
279 // OR
280 // 2) least squares linear fit to the at-most 24 nearest neighbors, with the function
281 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
282 // where the data is assume to live, are the same as cell centers. Note that calc_slopes_eb_grown
283 // defines the matrix A and passes this A to amrex_calc_slopes_eb_given_A_grown.
284 //
285 // This routine assumes that there are no relevant hoextrap/extdir domain boundary conditions for this cell --
286 // it does not test for them so this should not be called if such boundaries might be present
287 //
290 amrex_calc_slopes_eb_grown (int i, int j, int k, int n, int nx, int ny,
295  int max_order) noexcept
296 {
297  AMREX_ASSERT(nx == 1 || nx == 2);
298  AMREX_ASSERT(ny == 1 || ny == 2);
299 
300  constexpr int dim_a = 25;
301  amrex::Real A[dim_a][AMREX_SPACEDIM];
302 
303  int kk = 0;
304  {
305  // Make sure to zero all the entries in A (since the loop below may not cover all 25)
306  int lc=0;
307  for(int jj(-2); jj<=2; jj++){
308  for(int ii(-2); ii<=2; ii++){
309  A[lc][0] = amrex::Real(0.0);
310  A[lc][1] = amrex::Real(0.0);
311  lc++;
312  }
313  }
314 
315  lc=0;
316  for(int jj(-ny); jj<=ny; jj++){
317  for(int ii(-nx); ii<=nx; ii++){
318 
319  if(!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0))
320  {
321  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) - ccent(i,j,k,0);
322  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) - ccent(i,j,k,1);
323  }
324  lc++;
325  }
326  }
327  }
328 
329  const auto& eb_slopes = amrex_calc_slopes_eb_given_A_grown (i,j,k,n,nx,ny,A,state,flag);
330 
331  amrex::Real xslope = eb_slopes[0];
332  amrex::Real yslope = eb_slopes[1];
333 
334  // This will over-write the values of xslope and yslope if appropriate
335  amrex_overwrite_with_regular_slopes(i,j,k,n,xslope,yslope,state,vfrac,max_order);
336 
337  return {xslope,yslope};
338 }
339 
340 //
341 // amrex_overwrite_with_regular_slopes_extdir calculates the slope in each coordinate direction
342 // with a standard non-EB slope calculation (that depends on max_order)
343 //
345 void
347  amrex::Real& xslope, amrex::Real& yslope,
350  bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y,
351  int domlo_x, int domlo_y, int domhi_x, int domhi_y,
352  int max_order) noexcept
353 {
354  //
355  // Correct only those cells which are affected by extdir but not by EB:
356  // 2) If all the cells are regular we use the "regular slope" in the extdir direction
357  //
358 
359  // We have enough cells in the x-direction to do 4th order slopes
360  // centered on (i,j,k) with all values at cell centers
361  if (max_order == 4 &&
362  vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i-2,j,k) == 1. &&
363  vfrac(i+1,j,k) == 1. && vfrac(i+2,j,k) == 1.)
364  {
365  int order = 4;
366  xslope = amrex_calc_xslope_extdir(i,j,k,n,order,state,edlo_x,edhi_x,domlo_x,domhi_x);
367  }
368  // We have enough cells in the x-direction to do 2nd order slopes
369  // centered on (i,j,k) with all values at cell centers
370  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.)
371  {
372  int order = 2;
373  xslope = amrex_calc_xslope_extdir(i,j,k,n,order,state,edlo_x,edhi_x,domlo_x,domhi_x);
374  } else if (max_order == 0) {
375  xslope = 0.;
376  }
377 
378  // We have enough cells in the y-direction to do 4th order slopes
379  // centered on (i,j,k) with all values at cell centers
380  if (max_order == 4 &&
381  vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j-2,k) == 1. &&
382  vfrac(i,j+1,k) == 1. && vfrac(i,j+2,k) == 1.)
383  {
384  int order = 4;
385  yslope = amrex_calc_yslope_extdir(i,j,k,n,order,state,edlo_y,edhi_y,domlo_y,domhi_y);
386  }
387  // We have enough cells in the y-direction to do 2nd order slopes
388  // centered on (i,j,k) with all values at cell centers
389  else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.)
390  {
391  int order = 2;
392  yslope = amrex_calc_yslope_extdir(i,j,k,n,order,state,edlo_y,edhi_y,domlo_y,domhi_y);
393  } else if (max_order == 0) {
394  yslope = 0.;
395  }
396 }
397 
398 // amrex_calc_slopes_extdir_eb calculates the slope in each coordinate direction using a
399 // 1) standard limited slope if all three cells in the stencil are regular cells
400 // (this stencil sees the extdir/hoextrap boundary condition if there is one)
401 // OR
402 // 2) least squares linear fit to the at-most 8 nearest neighbors, with the function
403 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
404 // where the data is assume to live, are the same as cell centers.
405 //
408 amrex_calc_slopes_extdir_eb (int i, int j, int k, int n,
415  bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y,
416  int domlo_x, int domlo_y, int domhi_x, int domhi_y,
417  int max_order) noexcept
418 {
419  constexpr int dim_a = 9;
420 
421  auto xslope = amrex::Real(0.0);
422  auto yslope = amrex::Real(0.0);
423 
424  // First get EB-aware slope that doesn't know about extdir
425  bool needs_bdry_stencil = (edlo_x && i <= domlo_x) || (edhi_x && i >= domhi_x) ||
426  (edlo_y && j <= domlo_y) || (edhi_y && j >= domhi_y);
427 
428  //
429  // This call does not have any knowledge of extdir / hoextrap boundary conditions
430  //
431  if (!needs_bdry_stencil)
432  {
433  // This returns slopes calculated with the regular 1-d approach if all cells in the stencil
434  // are regular. If not, it uses the EB-aware least squares approach to fit a linear profile
435  // using the neighboring un-covered cells.
436  const auto& slopes = amrex_calc_slopes_eb (i,j,k,n,state,ccent,vfrac,flag,max_order);
437  return slopes;
438 
439  } else {
440 
441  amrex::Real A[dim_a][AMREX_SPACEDIM];
442 
443  int lc=0;
444  int kk = 0;
445 
446  for(int jj(-1); jj<=1; jj++) {
447  for(int ii(-1); ii<=1; ii++)
448  {
449  if( flag(i,j,k).isConnected(ii,jj,kk) &&
450  ! (ii==0 && jj==0 && kk==0))
451  {
452  bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1);
453  bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1);
454 
455  bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1);
456  bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1);
457 
458  bool klo_test = false;
459  bool khi_test = false;
460 
461  // These are the default values if no physical boundary
462  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0);
463  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1);
464  // Do corrections for entire x-face
465  if (ilo_test)
466  {
467  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
468  {
469  A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0);
470  }
471  A[lc][0] = amrex::Real(-0.5) ;
472  } else if (ihi_test) {
473 
474  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
475  {
476  A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0);
477  }
478  A[lc][0] = amrex::Real(0.5) ;
479  }
480 
481  // Do corrections for entire y-face
482  if (jlo_test)
483  {
484  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
485  {
486  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0);
487  }
488  A[lc][1] = amrex::Real(-0.5) ;
489 
490  } else if (jhi_test) {
491 
492  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
493  {
494  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0);
495  }
496  A[lc][1] = amrex::Real(0.5) ;
497  }
498 
499  A[lc][0] -= ccent(i,j,k,0);
500  A[lc][1] -= ccent(i,j,k,1);
501 
502  } else {
503  A[lc][0] = amrex::Real(0.0);
504  A[lc][1] = amrex::Real(0.0);
505  }
506  lc++;
507  }} // i,j
508 
509  const auto& slopes = amrex_calc_slopes_eb_given_A (i,j,k,n,A,state,flag);
510  xslope = slopes[0];
511  yslope = slopes[1];
512 
513  // This will over-write the values of xslope and yslope if appropriate
514  amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,state,vfrac,
515  edlo_x,edlo_y,edhi_x,edhi_y,
516  domlo_x,domlo_y,domhi_x,domhi_y,max_order);
517 
518  // Zero out slopes outside of an extdir (or hoextrap) boundary
519  // TODO: is this the right thing to do at a HOEXTRAP boundary??
520  if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) ||
521  (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) )
522  {
523  xslope = 0.;
524  yslope = 0.;
525  }
526  return {xslope,yslope};
527  }
528 }
529 
530 // amrex_calc_slopes_extdir_eb_grown calculates the slope in each coordinate direction using a
531 // 1) standard limited slope if all three cells in the stencil are regular cells
532 // (this stencil sees the extdir/hoextrap boundary condition if there is one)
533 // OR
534 // 2) least squares linear fit to the at-most 24 nearest neighbors, with the function
535 // going through the centroid of cell(i,j,k). This does not assume that the cell centroids,
536 // where the data is assume to live, are the same as cell centers.
537 //
540 amrex_calc_slopes_extdir_eb_grown (int i, int j, int k, int n, int nx, int ny,
547  bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y,
548  int domlo_x, int domlo_y, int domhi_x, int domhi_y,
549  int max_order) noexcept
550 {
551  constexpr int dim_a = 25;
552 
553  auto xslope = amrex::Real(0.0);
554  auto yslope = amrex::Real(0.0);
555 
556  // First get EB-aware slope that doesn't know about extdir
557  bool needs_bdry_stencil = (edlo_x && i <= domlo_x) || (edhi_x && i >= domhi_x) ||
558  (edlo_y && j <= domlo_y) || (edhi_y && j >= domhi_y);
559 
560  //
561  // This call does not have any knowledge of extdir / hoextrap boundary conditions
562  //
563  if (!needs_bdry_stencil)
564  {
565  // This returns slopes calculated with the regular 1-d approach if all cells in the stencil
566  // are regular. If not, it uses the EB-aware least squares approach to fit a linear profile
567  // using the neighboring un-covered cells.
568  const auto& slopes = amrex_calc_slopes_eb_grown (i,j,k,n,nx,ny,state,ccent,vfrac,flag,max_order);
569  return slopes;
570 
571  } else {
572 
573  amrex::Real A[dim_a][AMREX_SPACEDIM];
574 
575  int lc=0;
576  int kk = 0;
577 
578  for(int jj(-ny); jj<=ny; jj++) {
579  for(int ii(-nx); ii<=nx; ii++)
580  {
581  if ( !flag(i+ii,j+jj,k).isCovered() && !(ii==0 && jj==0 && kk==0) )
582  {
583  bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1);
584  bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1);
585 
586  bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1);
587  bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1);
588 
589  bool klo_test = false;
590  bool khi_test = false;
591 
592  // These are the default values if no physical boundary
593  A[lc][0] = amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0);
594  A[lc][1] = amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1);
595  // Do corrections for entire x-face
596  if (ilo_test)
597  {
598  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
599  {
600  A[lc][1] = amrex::Real(jj) + fcx(i ,j+jj,k+kk,0);
601  }
602  A[lc][0] = amrex::Real(-0.5) ;
603  } else if (ihi_test) {
604 
605  if (!jlo_test && !jhi_test && !klo_test && !khi_test)
606  {
607  A[lc][1] = amrex::Real(jj) + fcx(i+ii,j+jj,k+kk,0);
608  }
609  A[lc][0] = amrex::Real(0.5) ;
610  }
611 
612  // Do corrections for entire y-face
613  if (jlo_test)
614  {
615  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
616  {
617  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j ,k+kk,0);
618  }
619  A[lc][1] = amrex::Real(-0.5) ;
620 
621  } else if (jhi_test) {
622 
623  if (!ilo_test && !ihi_test && !klo_test && !khi_test)
624  {
625  A[lc][0] = amrex::Real(ii) + fcy(i+ii,j+jj,k+kk,0);
626  }
627  A[lc][1] = amrex::Real(0.5) ;
628  }
629 
630  A[lc][0] -= ccent(i,j,k,0);
631  A[lc][1] -= ccent(i,j,k,1);
632 
633  } else {
634  A[lc][0] = amrex::Real(0.0);
635  A[lc][1] = amrex::Real(0.0);
636  }
637  lc++;
638  }} // i,j
639 
640  const auto& slopes = amrex_calc_slopes_eb_given_A_grown (i,j,k,n,nx,ny,A,state,flag);
641  xslope = slopes[0];
642  yslope = slopes[1];
643 
644  // This will over-write the values of xslope and yslope if appropriate
645  amrex_overwrite_with_regular_slopes_extdir(i,j,k,n,xslope,yslope,state,vfrac,
646  edlo_x,edlo_y,edhi_x,edhi_y,
647  domlo_x,domlo_y,domhi_x,domhi_y,max_order);
648 
649  // Zero out slopes outside of an extdir (or hoextrap) boundary
650  // TODO: is this the right thing to do at a HOEXTRAP boundary??
651  if ( (edlo_x && i < domlo_x) || (edhi_x && i > domhi_x) ||
652  (edlo_y && j < domlo_y) || (edhi_y && j > domhi_y) )
653  {
654  xslope = 0.;
655  yslope = 0.;
656  }
657  return {xslope,yslope};
658  }
659 }
660 
662 amrex::Real
663 amrex_calc_alpha_stencil(amrex::Real q_hat, amrex::Real q_max,
664  amrex::Real q_min, amrex::Real state, amrex::Real alpha) noexcept
665 {
666  using namespace amrex::literals;
667 
668  auto alpha_temp = 0.0_rt;
669  auto small = 1.0e-13_rt;
670 
671  if ((q_hat-state) > small) {
672  alpha_temp = amrex::min(1.0_rt,(q_max-state)/(q_hat-state));
673  } else if ((q_hat-state) < -small) {
674  alpha_temp = amrex::min(1.0_rt,(q_min-state)/(q_hat-state));
675  } else {
676  alpha_temp = 1.0_rt;
677  }
678  return amrex::min(alpha, alpha_temp);
679 }
680 
683 amrex_calc_alpha_limiter(int i, int j, int k, int n,
689  amrex::Array4<amrex::Real const> const& ccent) noexcept
690 {
691  amrex::Real alpha = 2.0;
692 
693  int cuts_x = 0;
694  int cuts_y = 0;
695 
696  // Compute how many cut or regular faces do we have in 3x3 block
697  int kk = 0;
698  for(int jj(-1); jj<=1; jj++){
699  for(int ii(-1); ii<=1; ii++){
700  if( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0))
701  {
702  if ((ii==-1 || ii==1) && jj==0) { cuts_x++; }
703  if ((jj==-1 || jj==1) && ii==0) { cuts_y++; }
704  }
705  }
706  }
707 
708  amrex::Real xc = ccent(i,j,k,0); // centroid of cell (i,j,k)
709  amrex::Real yc = ccent(i,j,k,1);
710 
711  //Reconstruct values at the face centroids and compute the limiter
712  if(flag(i,j,k).isConnected(0,1,0)) {
713  amrex::Real xf = fcy(i,j+1,k,0); // local (x,z) of centroid of y-face we are extrapolating to
714 
715  amrex::Real delta_x = xf - xc;
716  amrex::Real delta_y = amrex::Real(0.5) - yc;
717 
718  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + delta_y * slopes[1];
719 
720  amrex::Real q_min = state(i,j,k,n);
721  amrex::Real q_max = state(i,j,k,n);
722 
723  // Compute max and min values in a 3x2 stencil
724  for(int jj(0); jj<=1; jj++){
725  for(int ii(-1); ii<=1; ii++){
726  if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) {
727  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
728  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
729  }
730  }
731  }
732 
733  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
734  }
735  if (flag(i,j,k).isConnected(0,-1,0)){
736  amrex::Real xf = fcy(i,j,k,0); // local (x,z) of centroid of y-face we are extrapolating to
737 
738  amrex::Real delta_x = xf - xc;
739  amrex::Real delta_y = amrex::Real(0.5) + yc;
740 
741  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] - delta_y * slopes[1];
742 
743  amrex::Real q_min = state(i,j,k,n);
744  amrex::Real q_max = state(i,j,k,n);
745 
746  // Compute max and min values in a 3x2 stencil
747  for(int jj(-1); jj<=0; jj++){
748  for(int ii(-1); ii<=1; ii++){
749  if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) {
750  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
751  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
752  }
753  }
754  }
755 
756  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
757  }
758  if(flag(i,j,k).isConnected(1,0,0)) {
759  amrex::Real yf = fcx(i+1,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to
760 
761  amrex::Real delta_x = amrex::Real(0.5) - xc;
762  amrex::Real delta_y = yf - yc;
763 
764  amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0] + delta_y * slopes[1];
765 
766  amrex::Real q_min = state(i,j,k,n);
767  amrex::Real q_max = state(i,j,k,n);
768 
769  for(int jj(-1); jj<=1; jj++){
770  for(int ii(0); ii<=1; ii++){
771  if ( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0) ) {
772  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
773  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
774  }
775  }
776  }
777 
778  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
779  }
780  if(flag(i,j,k).isConnected(-1,0,0)) {
781  amrex::Real yf = fcx(i,j,k,0); // local (y,z) of centroid of x-face we are extrapolating to
782 
783  amrex::Real delta_x = amrex::Real(0.5) + xc;
784  amrex::Real delta_y = yf - yc;
785 
786  amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0] + delta_y * slopes[1];
787 
788  amrex::Real q_min = state(i,j,k,n);
789  amrex::Real q_max = state(i,j,k,n);
790 
791  for(int jj(-1); jj<=1; jj++){
792  for(int ii(-1); ii<=0; ii++){
793  if( flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0)) {
794  if (state(i+ii,j+jj,k+kk,n) > q_max) { q_max = state(i+ii,j+jj,k+kk,n); }
795  if (state(i+ii,j+jj,k+kk,n) < q_min) { q_min = state(i+ii,j+jj,k+kk,n); }
796  }
797  }
798  }
799  alpha = amrex_calc_alpha_stencil(q_hat, q_max, q_min, state(i,j,k,n), alpha);
800  }
801 
802  amrex::Real xalpha = alpha;
803  amrex::Real yalpha = alpha;
804 
805  //Zeroing out the slopes in the direction where a covered face exists.
806  if (cuts_x<2) { xalpha = 0; }
807  if (cuts_y<2) { yalpha = 0; }
808 
809  return {xalpha,yalpha};
810 }
811 
812 //amrex_lim_slopes_eb computes the slopes calling amrex_calc_slopes_eb, and then each slope component
813 //is multiplied by a limiter based on the work of Barth-Jespersen.
816 amrex_lim_slopes_eb (int i, int j, int k, int n,
823  int max_order) noexcept
824 {
825 
828 
829  slopes = amrex_calc_slopes_eb(i,j,k,n,state,ccent,vfrac,flag,max_order);
830 
831  alpha_lim = amrex_calc_alpha_limiter(i,j,k,n,state,flag,slopes,fcx,fcy,ccent);
832 
833  // Setting limiter to 1 for stencils that just consists of non-EB cells because
834  // amrex_calc_slopes_eb routine will call the slope routine for non-EB stencils that has already a limiter
835  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) {
836  alpha_lim[0] = 1.0;
837  }
838 
839  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) {
840  alpha_lim[1] = 1.0;
841  }
842 
843  return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1]};
844 }
845 
846 //amrex_lim_slopes_extdir_eb computes the slopes calling amrex_calc_slopes_extdir_eb, and then each slope component
847 //is multiplied by a limiter based on the work of Barth-Jespersen.
850 amrex_lim_slopes_extdir_eb (int i, int j, int k, int n,
857  bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y,
858  int domlo_x, int domlo_y, int domhi_x, int domhi_y,
859  int max_order) noexcept
860 {
861 
864 
865  slopes = amrex_calc_slopes_extdir_eb(i,j,k,n,state,ccent,vfrac,fcx,fcy,flag,
866  edlo_x,edlo_y,edhi_x,edhi_y,
867  domlo_x,domlo_y,domhi_x,domhi_y,max_order);
868  alpha_lim = amrex_calc_alpha_limiter(i,j,k,n,state,flag,slopes,fcx,fcy,ccent);
869 
870  // Setting limiter to 1 for stencils that just consists of non-EB cells because
871  // amrex_calc_slopes_extdir_eb routine will call the slope routine for non-EB stencils that has already a limiter
872  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) {
873  alpha_lim[0] = 1.0;
874  }
875 
876  if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) {
877  alpha_lim[1] = 1.0;
878  }
879 
880  return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1]};
881 }
882 
883 #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, int n, int nx, int ny, amrex::Real A[25][AMREX_SPACEDIM], amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::EBCellFlag const > const &flag) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:96
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::EBCellFlag const > const &flag, int max_order) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:816
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_2D_K.H:663
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 &ccent) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:683
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::EBCellFlag const > const &flag, bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y, int domlo_x, int domlo_y, int domhi_x, int domhi_y, int max_order) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:408
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > amrex_calc_slopes_eb_given_A(int i, int j, int, int n, amrex::Real A[9][AMREX_SPACEDIM], amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::EBCellFlag const > const &flag) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:20
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::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &vfrac, int max_order) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:167
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::EBCellFlag const > const &flag, bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y, int domlo_x, int domlo_y, int domhi_x, int domhi_y, int max_order) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:850
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, 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::EBCellFlag const > const &flag, bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y, int domlo_x, int domlo_y, int domhi_x, int domhi_y, int max_order) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:540
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, 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_2D_K.H:290
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_2D_K.H:230
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::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &vfrac, bool edlo_x, bool edlo_y, bool edhi_x, bool edhi_y, int domlo_x, int domlo_y, int domhi_x, int domhi_y, int max_order) noexcept
Definition: AMReX_EB_Slopes_2D_K.H:346
#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