Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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//
20amrex_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//
96amrex_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//
166void
167amrex_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//
230amrex_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//
290amrex_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//
345void
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//
408amrex_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//
540amrex_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
662amrex::Real
663amrex_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 sml = 1.0e-13_rt;
670
671 if ((q_hat-state) > sml) {
672 alpha_temp = amrex::min(1.0_rt,(q_max-state)/(q_hat-state));
673 } else if ((q_hat-state) < -sml) {
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
683amrex_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.
816amrex_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.
850amrex_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_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(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 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_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 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_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_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_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 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 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
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_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
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:21
Definition AMReX_Array4.H:61
Definition AMReX_Array.H:34