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