ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_EBSlopes.H File Reference
#include <AMReX_EBCellFlag.H>
#include <ERF_IndexDefines.H>
Include dependency graph for ERF_EBSlopes.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet ([[maybe_unused]] int igrid_query, [[maybe_unused]] int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::RealVect const &bcent_eb, amrex::Real const state_eb, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet_staggered (int igrid_query, int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::RealVect const &bcent_eb, amrex::Real const state_eb, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_staggered (int igrid_query, int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
 
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_staggered_upwind (int igrid_query, [[maybe_unused]] int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, const amrex::Array4< const amrex::Real > &vel_arr, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
 

Function Documentation

◆ erf_calc_slopes_eb_Dirichlet()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> erf_calc_slopes_eb_Dirichlet ( [[maybe_unused] ] int  igrid_query,
[[maybe_unused] ] int  igrid_data,
amrex::Real  dx,
amrex::Real  dy,
amrex::Real  dz,
int  i,
int  j,
int  k,
amrex::RealVect const &  bcent_eb,
amrex::Real const  state_eb,
amrex::Array4< amrex::Real const > const &  state,
amrex::Array4< amrex::Real const > const &  ccent,
amrex::Array4< amrex::EBCellFlag const > const &  flag 
)

Function for computing the slopes with Dirichlet data This function works for only face-centered grids

Parameters
[in]igrid_queryindex of grids where the query point is located (0=CC, 1=U, 2=V, 3=W)
[in]igrid_dataindex of grids where data points are located (0=CC, 1=U, 2=V, 3=W)
[in]bcent_ebcentroid of the EB face
[in]state_ebstate at the EB face (Dirichlet data)
[in]statearray of state variable
[in]ccentarray of cell centroid coordinates
29 {
30  AMREX_ASSERT( igrid_query == igrid_data );
31 
32  // Fitting stencil
33 
34  amrex::RealVect bias{0.,0.,0.};
35 
36  constexpr int dim_a = 27;
37 
38  int ii_lo = -1;
39  int jj_lo = -1;
40  int kk_lo = -1;
41  int ii_hi = 1;
42  int jj_hi = 1;
43  int kk_hi = 1;
44 
45  amrex::Real A[dim_a][AMREX_SPACEDIM];
46 
47  // Array of distances from the query point
48  int ll=0;
49  for(int kk(kk_lo); kk<=kk_hi; kk++) {
50  for(int jj(jj_lo); jj<=jj_hi; jj++) {
51  for(int ii(ii_lo); ii<=ii_hi; ii++) {
52 
53  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
54  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
55  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
56  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
57  } else {
58  A[ll][0] = amrex::Real(0.0);
59  A[ll][1] = amrex::Real(0.0);
60  A[ll][2] = amrex::Real(0.0);
61  }
62  ll++;
63  }}}
64 
65  //
66  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
67  //
68  amrex::Real du[dim_a];
69 
70  ll=0;
71  for(int kk(kk_lo); kk<=kk_hi; kk++) {
72  for(int jj(jj_lo); jj<=jj_hi; jj++) {
73  for(int ii(ii_lo); ii<=ii_hi; ii++) {
74 
75  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
76  du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
77  } else {
78  du[ll] = amrex::Real(0.0);
79  }
80  ll++;
81  }}}
82 
83  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
84  amrex::Real Atb[AMREX_SPACEDIM];
85 
86  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
87  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
88  AtA[jj][ii] = amrex::Real(0.0);
89  }
90  Atb[jj] = amrex::Real(0.0);
91  }
92 
93  for(int lc(0); lc < dim_a; ++lc)
94  {
95  AtA[0][0] += A[lc][0]* A[lc][0];
96  AtA[0][1] += A[lc][0]* A[lc][1];
97  AtA[0][2] += A[lc][0]* A[lc][2];
98  AtA[1][1] += A[lc][1]* A[lc][1];
99  AtA[1][2] += A[lc][1]* A[lc][2];
100  AtA[2][2] += A[lc][2]* A[lc][2];
101 
102  Atb[0] += A[lc][0]*du[lc];
103  Atb[1] += A[lc][1]*du[lc];
104  Atb[2] += A[lc][2]*du[lc];
105  }
106 
107  // Fill in symmetric
108  AtA[1][0] = AtA[0][1];
109  AtA[2][0] = AtA[0][2];
110  AtA[2][1] = AtA[1][2];
111 
112  amrex::Real detAtA =
113  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
114  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
115  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
116 
117  amrex::Real detAtA_x =
118  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
119  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
120  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
121 
122  // Slope at centroid of (i,j,k)
123  amrex::Real xslope = detAtA_x / detAtA;
124 
125  amrex::Real detAtA_y =
126  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
127  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
128  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
129 
130  // Slope at centroid of (i,j,k)
131  amrex::Real yslope = detAtA_y / detAtA;
132 
133  amrex::Real detAtA_z =
134  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
135  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
136  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
137 
138  // Slope at centroid of (i,j,k)
139  amrex::Real zslope = detAtA_z / detAtA;
140 
141  return {xslope,yslope,zslope};
142 }
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:16

Referenced by DiffusionSrcForMom_EB().

Here is the caller graph for this function:

◆ erf_calc_slopes_eb_Dirichlet_staggered()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> erf_calc_slopes_eb_Dirichlet_staggered ( int  igrid_query,
int  igrid_data,
amrex::Real  dx,
amrex::Real  dy,
amrex::Real  dz,
int  i,
int  j,
int  k,
amrex::RealVect const &  bcent_eb,
amrex::Real const  state_eb,
amrex::Array4< amrex::Real const > const &  state,
amrex::Array4< amrex::Real const > const &  ccent,
amrex::Array4< amrex::EBCellFlag const > const &  flag 
)

Function for computing the slopes with Dirichlet data in staggered grids This function works for only face-centered grids

Parameters
[in]igrid_queryindex of grids where the query point is located (0=CC, 1=U, 2=V, 3=W)
[in]igrid_dataindex of grids where data points are located (0=CC, 1=U, 2=V, 3=W)
[in]bcent_ebcentroid of the EB face
[in]state_ebstate at the EB face (Dirichlet data)
[in]statearray of state variable
[in]ccentarray of cell centroid coordinates
166 {
167  // Fitting stencil
168 
169  constexpr int dim_a = 48;
170 
171  int ii_lo = -1;
172  int jj_lo = -1;
173  int kk_lo = -1;
174  int ii_hi = 1;
175  int jj_hi = 1;
176  int kk_hi = 1;
177 
178  if (igrid_query == Vars::xvel) {
179  ii_lo = -2;
180  } else if (igrid_query == Vars::yvel) {
181  jj_lo = -2;
182  } else if (igrid_query == Vars::zvel) {
183  kk_lo = -2;
184  }
185 
186  if (igrid_data == Vars::xvel) {
187  ii_hi = 2;
188  } else if (igrid_data == Vars::yvel) {
189  jj_hi = 2;
190  } else if (igrid_data == Vars::zvel) {
191  kk_hi = 2;
192  }
193 
194  // Set bias in the index space between the two staggered grids
195  amrex::RealVect bias{0.,0.,0.};
196  bias[igrid_query-1] = 0.5;
197  bias[igrid_data-1] = -0.5;
198 
199  amrex::Real A[dim_a][AMREX_SPACEDIM];
200 
201  // Array of distances from the query point
202  int ll=0;
203  for(int kk(kk_lo); kk<=kk_hi; kk++) {
204  for(int jj(jj_lo); jj<=jj_hi; jj++) {
205  for(int ii(ii_lo); ii<=ii_hi; ii++) {
206 
207  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
208  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
209  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
210  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
211  } else {
212  A[ll][0] = amrex::Real(0.0);
213  A[ll][1] = amrex::Real(0.0);
214  A[ll][2] = amrex::Real(0.0);
215  }
216  ll++;
217  }}}
218 
219  //
220  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
221  //
222  amrex::Real du[dim_a];
223 
224  ll=0;
225  for(int kk(kk_lo); kk<=kk_hi; kk++) {
226  for(int jj(jj_lo); jj<=jj_hi; jj++) {
227  for(int ii(ii_lo); ii<=ii_hi; ii++) {
228 
229  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
230  du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
231  } else {
232  du[ll] = amrex::Real(0.0);
233  }
234  ll++;
235  }}}
236 
237  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
238  amrex::Real Atb[AMREX_SPACEDIM];
239 
240  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
241  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
242  AtA[jj][ii] = amrex::Real(0.0);
243  }
244  Atb[jj] = amrex::Real(0.0);
245  }
246 
247  for(int lc(0); lc < dim_a; ++lc)
248  {
249  AtA[0][0] += A[lc][0]* A[lc][0];
250  AtA[0][1] += A[lc][0]* A[lc][1];
251  AtA[0][2] += A[lc][0]* A[lc][2];
252  AtA[1][1] += A[lc][1]* A[lc][1];
253  AtA[1][2] += A[lc][1]* A[lc][2];
254  AtA[2][2] += A[lc][2]* A[lc][2];
255 
256  Atb[0] += A[lc][0]*du[lc];
257  Atb[1] += A[lc][1]*du[lc];
258  Atb[2] += A[lc][2]*du[lc];
259  }
260 
261  // Fill in symmetric
262  AtA[1][0] = AtA[0][1];
263  AtA[2][0] = AtA[0][2];
264  AtA[2][1] = AtA[1][2];
265 
266  amrex::Real detAtA =
267  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
268  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
269  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
270 
271  amrex::Real detAtA_x =
272  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
273  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
274  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
275 
276  // Slope at centroid of (i,j,k)
277  amrex::Real xslope = detAtA_x / detAtA;
278 
279  amrex::Real detAtA_y =
280  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
281  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
282  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
283 
284  // Slope at centroid of (i,j,k)
285  amrex::Real yslope = detAtA_y / detAtA;
286 
287  amrex::Real detAtA_z =
288  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
289  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
290  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
291 
292  // Slope at centroid of (i,j,k)
293  amrex::Real zslope = detAtA_z / detAtA;
294 
295  return {xslope,yslope,zslope};
296 }
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:143
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142

Referenced by DiffusionSrcForMom_EB().

Here is the caller graph for this function:

◆ erf_calc_slopes_eb_staggered()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> erf_calc_slopes_eb_staggered ( int  igrid_query,
int  igrid_data,
amrex::Real  dx,
amrex::Real  dy,
amrex::Real  dz,
int  i,
int  j,
int  k,
amrex::Array4< amrex::Real const > const &  state,
amrex::Array4< amrex::Real const > const &  ccent,
amrex::Array4< amrex::EBCellFlag const > const &  flag 
)
307 {
308  AMREX_ASSERT((igrid_query == Vars::xvel || igrid_query == Vars::yvel || igrid_query == Vars::zvel) && igrid_data == Vars::cons);
309 
310  // Fitting stencil
311 
312  constexpr int dim_a = 36;
313 
314  int ii_lo = -1;
315  int jj_lo = -1;
316  int kk_lo = -1;
317  int ii_hi = 1;
318  int jj_hi = 1;
319  int kk_hi = 1;
320 
321  if (igrid_query == Vars::xvel) {
322  ii_lo = -2;
323  } else if (igrid_query == Vars::yvel) {
324  jj_lo = -2;
325  } else if (igrid_query == Vars::zvel) {
326  kk_lo = -2;
327  }
328 
329  if (igrid_data == Vars::xvel) {
330  ii_hi = 2;
331  } else if (igrid_data == Vars::yvel) {
332  jj_hi = 2;
333  } else if (igrid_data == Vars::zvel) {
334  kk_hi = 2;
335  }
336 
337  // Set bias in the index space between the two staggered grids
338  amrex::RealVect bias{0.,0.,0.};
339  bias[igrid_query-1] = 0.5;
340 
341  amrex::Real A[dim_a][AMREX_SPACEDIM];
342 
343  // Array of distances from the query point
344  int ll=0;
345  for(int kk(kk_lo); kk<=kk_hi; kk++) {
346  for(int jj(jj_lo); jj<=jj_hi; jj++) {
347  for(int ii(ii_lo); ii<=ii_hi; ii++) {
348 
349  // if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
350  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
351  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - ccent(i,j,k,0) ) * dx;
352  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - ccent(i,j,k,1) ) * dy;
353  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - ccent(i,j,k,2) ) * dz;
354  } else {
355  A[ll][0] = amrex::Real(0.0);
356  A[ll][1] = amrex::Real(0.0);
357  A[ll][2] = amrex::Real(0.0);
358  }
359  ll++;
360  }}}
361 
362  //
363  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
364  //
365  amrex::Real du[dim_a];
366 
367  ll=0;
368  for(int kk(kk_lo); kk<=kk_hi; kk++) {
369  for(int jj(jj_lo); jj<=jj_hi; jj++) {
370  for(int ii(ii_lo); ii<=ii_hi; ii++) {
371 
372  // if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
373  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
374  du[ll] = state(i+ii,j+jj,k+kk) - state(i,j,k);
375  } else {
376  du[ll] = amrex::Real(0.0);
377  }
378  ll++;
379  }}}
380 
381  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
382  amrex::Real Atb[AMREX_SPACEDIM];
383 
384  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
385  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
386  AtA[jj][ii] = amrex::Real(0.0);
387  }
388  Atb[jj] = amrex::Real(0.0);
389  }
390 
391  for(int lc(0); lc < dim_a; ++lc)
392  {
393  AtA[0][0] += A[lc][0]* A[lc][0];
394  AtA[0][1] += A[lc][0]* A[lc][1];
395  AtA[0][2] += A[lc][0]* A[lc][2];
396  AtA[1][1] += A[lc][1]* A[lc][1];
397  AtA[1][2] += A[lc][1]* A[lc][2];
398  AtA[2][2] += A[lc][2]* A[lc][2];
399 
400  Atb[0] += A[lc][0]*du[lc];
401  Atb[1] += A[lc][1]*du[lc];
402  Atb[2] += A[lc][2]*du[lc];
403  }
404 
405  // Fill in symmetric
406  AtA[1][0] = AtA[0][1];
407  AtA[2][0] = AtA[0][2];
408  AtA[2][1] = AtA[1][2];
409 
410  amrex::Real detAtA =
411  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
412  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
413  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
414 
415  amrex::Real detAtA_x =
416  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
417  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
418  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
419 
420  // Slope at centroid of (i,j,k)
421  amrex::Real xslope = detAtA_x / detAtA;
422 
423  amrex::Real detAtA_y =
424  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
425  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
426  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
427 
428  // Slope at centroid of (i,j,k)
429  amrex::Real yslope = detAtA_y / detAtA;
430 
431  amrex::Real detAtA_z =
432  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
433  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
434  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
435 
436  // Slope at centroid of (i,j,k)
437  amrex::Real zslope = detAtA_z / detAtA;
438 
439  return {xslope,yslope,zslope};
440 }
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140

Referenced by compute_gradp().

Here is the caller graph for this function:

◆ erf_calc_slopes_eb_staggered_upwind()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> erf_calc_slopes_eb_staggered_upwind ( int  igrid_query,
[[maybe_unused] ] int  igrid_data,
amrex::Real  dx,
amrex::Real  dy,
amrex::Real  dz,
int  i,
int  j,
int  k,
const amrex::Array4< const amrex::Real > &  vel_arr,
amrex::Array4< amrex::Real const > const &  state,
amrex::Array4< amrex::Real const > const &  ccent,
amrex::Array4< amrex::EBCellFlag const > const &  flag 
)
452 {
453  AMREX_ASSERT((igrid_query == Vars::xvel || igrid_query == Vars::yvel || igrid_query == Vars::zvel) && igrid_data == Vars::cons);
454 
455  // Fitting stencil
456 
457  int ii_lo = -1;
458  int jj_lo = -1;
459  int kk_lo = -1;
460  int ii_hi = 1;
461  int jj_hi = 1;
462  int kk_hi = 1;
463 
464  if (igrid_query == Vars::xvel && vel_arr(i,j,k)>0.0) {
465  ii_lo = -2;
466  ii_hi = 0;
467  }
468  if (igrid_query == Vars::yvel && vel_arr(i,j,k)>0.0) {
469  jj_lo = -2;
470  jj_hi = 0;
471  }
472  if (igrid_query == Vars::zvel && vel_arr(i,j,k)>0.0) {
473  kk_lo = -2;
474  kk_hi = 0;
475  }
476 
477  // Set bias in the index space between the two staggered grids
478  amrex::RealVect bias{0.,0.,0.};
479  bias[igrid_query-1] = 0.5;
480 
481  constexpr int dim_a = 27;
482 
483  amrex::Real A[dim_a][AMREX_SPACEDIM];
484 
485  // Array of distances from the query point
486  int ll=0;
487  for(int kk(kk_lo); kk<=kk_hi; kk++) {
488  for(int jj(jj_lo); jj<=jj_hi; jj++) {
489  for(int ii(ii_lo); ii<=ii_hi; ii++) {
490 
491  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
492  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - ccent(i,j,k,0) ) * dx;
493  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - ccent(i,j,k,1) ) * dy;
494  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - ccent(i,j,k,2) ) * dz;
495  } else {
496  A[ll][0] = amrex::Real(0.0);
497  A[ll][1] = amrex::Real(0.0);
498  A[ll][2] = amrex::Real(0.0);
499  }
500  ll++;
501  }}}
502 
503  //
504  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
505  //
506  amrex::Real du[dim_a];
507 
508  ll=0;
509  for(int kk(kk_lo); kk<=kk_hi; kk++) {
510  for(int jj(jj_lo); jj<=jj_hi; jj++) {
511  for(int ii(ii_lo); ii<=ii_hi; ii++) {
512 
513  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
514  du[ll] = state(i+ii,j+jj,k+kk) - state(i,j,k);
515  } else {
516  du[ll] = amrex::Real(0.0);
517  }
518  ll++;
519  }}}
520 
521  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
522  amrex::Real Atb[AMREX_SPACEDIM];
523 
524  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
525  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
526  AtA[jj][ii] = amrex::Real(0.0);
527  }
528  Atb[jj] = amrex::Real(0.0);
529  }
530 
531  for(int lc(0); lc < dim_a; ++lc)
532  {
533  AtA[0][0] += A[lc][0]* A[lc][0];
534  AtA[0][1] += A[lc][0]* A[lc][1];
535  AtA[0][2] += A[lc][0]* A[lc][2];
536  AtA[1][1] += A[lc][1]* A[lc][1];
537  AtA[1][2] += A[lc][1]* A[lc][2];
538  AtA[2][2] += A[lc][2]* A[lc][2];
539 
540  Atb[0] += A[lc][0]*du[lc];
541  Atb[1] += A[lc][1]*du[lc];
542  Atb[2] += A[lc][2]*du[lc];
543  }
544 
545  // Fill in symmetric
546  AtA[1][0] = AtA[0][1];
547  AtA[2][0] = AtA[0][2];
548  AtA[2][1] = AtA[1][2];
549 
550  amrex::Real detAtA =
551  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
552  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
553  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
554 
555  amrex::Real detAtA_x =
556  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
557  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
558  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
559 
560  // Slope at centroid of (i,j,k)
561  amrex::Real xslope = detAtA_x / detAtA;
562 
563  amrex::Real detAtA_y =
564  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
565  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
566  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
567 
568  // Slope at centroid of (i,j,k)
569  amrex::Real yslope = detAtA_y / detAtA;
570 
571  amrex::Real detAtA_z =
572  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
573  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
574  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
575 
576  // Slope at centroid of (i,j,k)
577  amrex::Real zslope = detAtA_z / detAtA;
578 
579  return {xslope,yslope,zslope};
580 }