1 #ifndef AMREX_HYDRO_EB_SLOPES_3D_K_H_
2 #define AMREX_HYDRO_EB_SLOPES_3D_K_H_
4 #include <AMReX_Config.H>
20 amrex::Real A[27][AMREX_SPACEDIM],
24 constexpr
int dim_a = 27;
25 amrex::Real du[dim_a];
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))
33 du[ll] = state(i+ii,j+jj,k+kk,n) - state(i,j,k,n);
35 du[ll] = amrex::Real(0.0);
40 amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
41 amrex::Real Atb[AMREX_SPACEDIM];
43 for(
int jj(0); jj<AMREX_SPACEDIM; ++jj){
44 for(
int ii(0); ii<AMREX_SPACEDIM; ++ii){
45 AtA[jj][ii] = amrex::Real(0.0);
47 Atb[jj] = amrex::Real(0.0);
50 for(
int lc(0); lc < 27; ++lc)
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];
59 Atb[0] += A[lc][0]*du[lc];
60 Atb[1] += A[lc][1]*du[lc];
61 Atb[2] += A[lc][2]*du[lc];
65 AtA[1][0] = AtA[0][1];
66 AtA[2][0] = AtA[0][2];
67 AtA[2][1] = AtA[1][2];
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]);
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] );
80 amrex::Real xs = detAtA_x / detAtA;
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]);
88 amrex::Real ys = detAtA_y / detAtA;
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]);
96 amrex::Real zs = detAtA_z / detAtA;
111 amrex::Real A[125][AMREX_SPACEDIM],
119 amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
120 amrex::Real Atb[AMREX_SPACEDIM];
122 for(
int jj(0); jj<AMREX_SPACEDIM; ++jj){
123 for(
int ii(0); ii<AMREX_SPACEDIM; ++ii){
124 AtA[jj][ii] = amrex::Real(0.0);
126 Atb[jj] = amrex::Real(0.0);
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];
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);
148 Atb[0] += A[lc][0]*du;
149 Atb[1] += A[lc][1]*du;
150 Atb[2] += A[lc][2]*du;
156 AtA[1][0] = AtA[0][1];
157 AtA[2][0] = AtA[0][2];
158 AtA[2][1] = AtA[1][2];
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]);
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] );
171 amrex::Real xs = detAtA_x / detAtA;
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]);
179 amrex::Real ys = detAtA_y / detAtA;
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]);
187 amrex::Real zs = detAtA_z / detAtA;
199 amrex::Real& xslope, amrex::Real& yslope,
203 int max_order) noexcept
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.)
221 else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.)
225 }
else if (max_order == 0) {
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.)
240 else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.)
244 }
else if (max_order == 0) {
248 #if (AMREX_SPACEDIM == 3)
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.)
256 zslope = amrex_calc_zslope(i,j,k,n,order,state);
260 else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.)
263 zslope = amrex_calc_zslope(i,j,k,n,order,state);
264 }
else if (max_order == 0) {
287 int max_order) noexcept
289 constexpr
int dim_a = 27;
290 amrex::Real A[dim_a][AMREX_SPACEDIM];
293 for(
int kk(-1); kk<=1; kk++) {
294 for(
int jj(-1); jj<=1; jj++) {
295 for(
int ii(-1); ii<=1; ii++)
297 if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0))
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);
303 A[lc][0] = amrex::Real(0.0);
304 A[lc][1] = amrex::Real(0.0);
305 A[lc][2] = amrex::Real(0.0);
314 amrex::Real xslope = slopes[0];
315 amrex::Real yslope = slopes[1];
316 amrex::Real zslope = slopes[2];
321 return {xslope,yslope,zslope};
341 int max_order) noexcept
347 constexpr
int dim_a = 125;
348 amrex::Real A[dim_a][AMREX_SPACEDIM];
352 for(
int kk(-2); kk<=2; kk++) {
353 for(
int jj(-2); jj<=2; jj++) {
354 for(
int ii(-2); ii<=2; ii++)
356 A[lc][0] = amrex::Real(0.0);
357 A[lc][1] = amrex::Real(0.0);
358 A[lc][2] = amrex::Real(0.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++)
367 if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0))
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);
379 amrex::Real xslope = slopes[0];
380 amrex::Real yslope = slopes[1];
381 amrex::Real zslope = slopes[2];
386 return {xslope,yslope,zslope};
396 amrex::Real& xslope, amrex::Real& yslope,
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
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.)
421 else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.)
425 }
else if (max_order == 0) {
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.)
440 else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.)
444 }
else if (max_order == 0) {
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.)
455 zslope = amrex_calc_zslope_extdir(i,j,k,n,order,state,edlo_z,edhi_z,domlo_z,domhi_z);
459 else if (max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.)
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) {
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
495 constexpr
int dim_a = 27;
497 auto xslope = amrex::Real(0.0);
498 auto yslope = amrex::Real(0.0);
499 auto zslope = amrex::Real(0.0);
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) ;
509 if (!needs_bdry_stencil)
519 amrex::Real A[dim_a][AMREX_SPACEDIM];
522 for(
int kk(-1); kk<=1; kk++) {
523 for(
int jj(-1); jj<=1; jj++) {
524 for(
int ii(-1); ii<=1; ii++)
526 if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0))
528 bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1);
529 bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1);
531 bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1);
532 bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1);
534 bool klo_test = ( edlo_z && (k == domlo_z) && kk == -1);
535 bool khi_test = ( edhi_z && (k == domhi_z) && kk == 1);
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);
544 if (!jlo_test && !jhi_test && !klo_test && !khi_test)
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);
549 A[lc][0] = -amrex::Real(0.5);
550 }
else if (ihi_test) {
552 if (!jlo_test && !jhi_test && !klo_test && !khi_test)
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);
557 A[lc][0] = amrex::Real(0.5);
563 if (!ilo_test && !ihi_test && !klo_test && !khi_test)
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);
568 A[lc][1] = -amrex::Real(0.5);
570 }
else if (jhi_test) {
572 if (!ilo_test && !ihi_test && !klo_test && !khi_test)
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);
577 A[lc][1] = amrex::Real(0.5);
583 if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
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);
588 A[lc][2] = -amrex::Real(0.5);
590 }
else if (khi_test) {
591 if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
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);
596 A[lc][2] = amrex::Real(0.5);
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);
605 A[lc][0] = amrex::Real(0.0);
606 A[lc][1] = amrex::Real(0.0);
607 A[lc][2] = amrex::Real(0.0);
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,
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) )
631 xslope = 0.; yslope = 0.; zslope = 0.;
634 return {xslope,yslope,zslope};
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
662 constexpr
int dim_a = 125;
664 auto xslope = amrex::Real(0.0);
665 auto yslope = amrex::Real(0.0);
666 auto zslope = amrex::Real(0.0);
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) ;
676 if (!needs_bdry_stencil)
681 const auto& slopes =
amrex_calc_slopes_eb_grown (i,j,k,n,nx,ny,nz,state,ccent,vfrac,flag,max_order);
686 amrex::Real A[dim_a][AMREX_SPACEDIM];
689 for(
int kk(-nz); kk<=nz; kk++) {
690 for(
int jj(-ny); jj<=ny; jj++) {
691 for(
int ii(-nx); ii<=nx; ii++)
693 if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0))
695 bool ilo_test = ( edlo_x && (i == domlo_x) && ii == -1);
696 bool ihi_test = ( edhi_x && (i == domhi_x) && ii == 1);
698 bool jlo_test = ( edlo_y && (j == domlo_y) && jj == -1);
699 bool jhi_test = ( edhi_y && (j == domhi_y) && jj == 1);
701 bool klo_test = ( edlo_z && (k == domlo_z) && kk == -1);
702 bool khi_test = ( edhi_z && (k == domhi_z) && kk == 1);
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);
711 if (!jlo_test && !jhi_test && !klo_test && !khi_test)
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);
716 A[lc][0] = -amrex::Real(0.5);
717 }
else if (ihi_test) {
719 if (!jlo_test && !jhi_test && !klo_test && !khi_test)
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);
724 A[lc][0] = amrex::Real(0.5);
730 if (!ilo_test && !ihi_test && !klo_test && !khi_test)
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);
735 A[lc][1] = -amrex::Real(0.5);
737 }
else if (jhi_test) {
739 if (!ilo_test && !ihi_test && !klo_test && !khi_test)
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);
744 A[lc][1] = amrex::Real(0.5);
750 if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
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);
755 A[lc][2] = -amrex::Real(0.5);
757 }
else if (khi_test) {
758 if (!ilo_test && !ihi_test && !jlo_test && !jhi_test)
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);
763 A[lc][2] = amrex::Real(0.5);
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);
771 A[lc][0] = amrex::Real(0.0);
772 A[lc][1] = amrex::Real(0.0);
773 A[lc][2] = amrex::Real(0.0);
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,
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) )
797 xslope = 0.; yslope = 0.; zslope = 0.;
800 return {xslope,yslope,zslope};
806 amrex::Real q_min, amrex::Real state, amrex::Real alpha) noexcept
808 using namespace amrex::literals;
810 auto alpha_temp = amrex::Real(0.0);
811 auto small = amrex::Real(1.0e-13);
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));
834 auto alpha = amrex::Real(2.0);
836 int cuts_x = 0;
int cuts_y = 0;
int cuts_z = 0;
839 for(
int kk(-1); kk<=1; kk++)
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))
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++; }
853 amrex::Real xc = ccent(i,j,k,0);
854 amrex::Real yc = ccent(i,j,k,1);
855 amrex::Real zc = ccent(i,j,k,2);
859 if (flag(i,j,k).isConnected(0,1,0)) {
860 amrex::Real xf = fcy(i,j+1,k,0);
861 amrex::Real zf = fcy(i,j+1,k,1);
863 amrex::Real delta_x = xf - xc;
864 amrex::Real delta_y = amrex::Real(0.5) - yc;
865 amrex::Real delta_z = zf - zc;
867 amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
868 + delta_y * slopes[1]
869 + delta_z * slopes[2];
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); }
887 if (flag(i,j,k).isConnected(0,-1,0)){
888 amrex::Real xf = fcy(i,j,k,0);
889 amrex::Real zf = fcy(i,j,k,1);
891 amrex::Real delta_x = xf - xc;
892 amrex::Real delta_y = amrex::Real(0.5) + yc;
893 amrex::Real delta_z = zf - zc;
895 amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
896 - delta_y * slopes[1]
897 + delta_z * slopes[2];
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); }
915 if (flag(i,j,k).isConnected(1,0,0)) {
916 amrex::Real yf = fcx(i+1,j,k,0);
917 amrex::Real zf = fcx(i+1,j,k,1);
919 amrex::Real delta_x = amrex::Real(0.5) - xc;
920 amrex::Real delta_y = yf - yc;
921 amrex::Real delta_z = zf - zc;
923 amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
924 + delta_y * slopes[1]
925 + delta_z * slopes[2];
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); }
943 if (flag(i,j,k).isConnected(-1,0,0)) {
944 amrex::Real yf = fcx(i,j,k,0);
945 amrex::Real zf = fcx(i,j,k,1);
947 amrex::Real delta_x = amrex::Real(0.5) + xc;
948 amrex::Real delta_y = yf - yc;
949 amrex::Real delta_z = zf - zc;
951 amrex::Real q_hat = state(i,j,k,n) - delta_x * slopes[0]
952 + delta_y * slopes[1]
953 + delta_z * slopes[2];
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); }
970 if (flag(i,j,k).isConnected(0,0,1)) {
971 amrex::Real xf = fcz(i,j,k+1,0);
972 amrex::Real yf = fcz(i,j,k+1,1);
974 amrex::Real delta_x = xf - xc;
975 amrex::Real delta_y = yf - yc;
976 amrex::Real delta_z = amrex::Real(0.5) - zc;
978 amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
979 + delta_y * slopes[1]
980 + delta_z * slopes[2];
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); }
997 if (flag(i,j,k).isConnected(0,0,-1)) {
998 amrex::Real xf = fcz(i,j,k,0);
999 amrex::Real yf = fcz(i,j,k,1);
1001 amrex::Real delta_x = xf - xc;
1002 amrex::Real delta_y = yf - yc;
1003 amrex::Real delta_z = amrex::Real(0.5) + zc;
1005 amrex::Real q_hat = state(i,j,k,n) + delta_x * slopes[0]
1006 + delta_y * slopes[1]
1007 - delta_z * slopes[2];
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); }
1025 amrex::Real xalpha = alpha;
1026 amrex::Real yalpha = alpha;
1027 amrex::Real zalpha = alpha;
1030 if (cuts_x<2) { xalpha = 0; }
1031 if (cuts_y<2) { yalpha = 0; }
1032 if (cuts_z<2) { zalpha = 0; }
1034 return {xalpha,yalpha,zalpha};
1049 int max_order) noexcept
1057 alpha_lim =
amrex_calc_alpha_limiter(i,j,k,n,state,flag,slopes,fcx,fcy,fcz,ccent);
1061 if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) {
1065 if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) {
1069 if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) {
1073 return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1],alpha_lim[2]*slopes[2]};
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
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);
1105 if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i-1,j,k) == 1. && vfrac(i+1,j,k) == 1.) {
1109 if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j-1,k) == 1. && vfrac(i,j+1,k) == 1.) {
1113 if ( max_order > 0 && vfrac(i,j,k) == 1. && vfrac(i,j,k-1) == 1. && vfrac(i,j,k+1) == 1.) {
1117 return {alpha_lim[0]*slopes[0],alpha_lim[1]*slopes[1],alpha_lim[2]*slopes[2]};
#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