ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_EBSlopes.H
Go to the documentation of this file.
1 #ifndef ERF_EBSLOPES_H_
2 #define ERF_EBSLOPES_H_
3 
4 #include <AMReX_EBCellFlag.H>
5 #include <ERF_IndexDefines.H>
6 
7 /**
8  * Function for computing the slopes with Dirichlet data
9  * This function works for only face-centered grids
10  *
11  * @param[in] igrid_query index of grids where the query point is located (0=CC, 1=U, 2=V, 3=W)
12  * @param[in] igrid_data index of grids where data points are located (0=CC, 1=U, 2=V, 3=W)
13  * @param[in] bcent_eb centroid of the EB face
14  * @param[in] state_eb state at the EB face (Dirichlet data)
15  * @param[in] state array of state variable
16  * @param[in] ccent array of cell centroid coordinates
17  */
18 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
19 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
20 erf_calc_slopes_eb_Dirichlet (
21  amrex::Real dx, amrex::Real dy, amrex::Real dz,
22  int i, int j, int k,
23  amrex::RealVect const& bcent_eb,
24  amrex::Real const state_eb,
25  amrex::Array4<amrex::Real const> const& state,
26  amrex::Array4<amrex::Real const> const& ccent,
27  amrex::Array4<amrex::EBCellFlag const> const& flag)
28 {
29  constexpr int dim_a = 27;
30 
31  int ii_lo = -1; int jj_lo = -1; int kk_lo = -1;
32  int ii_hi = 1; int jj_hi = 1; int kk_hi = 1;
33 
34  amrex::Real A[dim_a][AMREX_SPACEDIM];
35 
36  // Array of distances from the query point
37  int ll=0;
38  for(int kk(kk_lo); kk<=kk_hi; kk++) {
39  for(int jj(jj_lo); jj<=jj_hi; jj++) {
40  for(int ii(ii_lo); ii<=ii_hi; ii++) {
41 
42  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
43  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) - bcent_eb[0] ) * dx;
44  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) - bcent_eb[1] ) * dy;
45  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) - bcent_eb[2] ) * dz;
46  } else {
47  A[ll][0] = amrex::Real(0.0);
48  A[ll][1] = amrex::Real(0.0);
49  A[ll][2] = amrex::Real(0.0);
50  }
51  ll++;
52  }}}
53 
54  //
55  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
56  //
57  amrex::Real du[dim_a];
58 
59  ll=0;
60  for(int kk(kk_lo); kk<=kk_hi; kk++) {
61  for(int jj(jj_lo); jj<=jj_hi; jj++) {
62  for(int ii(ii_lo); ii<=ii_hi; ii++) {
63 
64  if (flag(i,j,k).isConnected(ii,jj,kk) && !(ii==0 && jj==0 && kk==0)) {
65  du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
66  } else {
67  du[ll] = amrex::Real(0.0);
68  }
69  ll++;
70  }}}
71 
72  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
73  amrex::Real Atb[AMREX_SPACEDIM];
74 
75  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
76  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
77  AtA[jj][ii] = amrex::Real(0.0);
78  }
79  Atb[jj] = amrex::Real(0.0);
80  }
81 
82  for(int lc(0); lc < dim_a; ++lc)
83  {
84  AtA[0][0] += A[lc][0]* A[lc][0];
85  AtA[0][1] += A[lc][0]* A[lc][1];
86  AtA[0][2] += A[lc][0]* A[lc][2];
87  AtA[1][1] += A[lc][1]* A[lc][1];
88  AtA[1][2] += A[lc][1]* A[lc][2];
89  AtA[2][2] += A[lc][2]* A[lc][2];
90 
91  Atb[0] += A[lc][0]*du[lc];
92  Atb[1] += A[lc][1]*du[lc];
93  Atb[2] += A[lc][2]*du[lc];
94  }
95 
96  // Fill in symmetric
97  AtA[1][0] = AtA[0][1];
98  AtA[2][0] = AtA[0][2];
99  AtA[2][1] = AtA[1][2];
100 
101  amrex::Real detAtA =
102  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
103  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
104  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
105 
106  amrex::Real detAtA_x =
107  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
108  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
109  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
110 
111  // Slope at centroid of (i,j,k)
112  amrex::Real xslope = detAtA_x / detAtA;
113 
114  amrex::Real detAtA_y =
115  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
116  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
117  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
118 
119  // Slope at centroid of (i,j,k)
120  amrex::Real yslope = detAtA_y / detAtA;
121 
122  amrex::Real detAtA_z =
123  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
124  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
125  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
126 
127  // Slope at centroid of (i,j,k)
128  amrex::Real zslope = detAtA_z / detAtA;
129 
130  return {xslope,yslope,zslope};
131 }
132 
133 /**
134  * Function for computing the slopes with Dirichlet data in staggered grids
135  * This function works for only face-centered grids
136  *
137  * @param[in] igrid_query index of grids where the query point is located (0=CC, 1=U, 2=V, 3=W)
138  * @param[in] igrid_data index of grids where data points are located (0=CC, 1=U, 2=V, 3=W)
139  * @param[in] bcent_eb centroid of the EB face
140  * @param[in] state_eb state at the EB face (Dirichlet data)
141  * @param[in] state array of state variable
142  * @param[in] ccent array of cell centroid coordinates
143  */
144 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
145 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
146 erf_calc_slopes_eb_Dirichlet_staggered (int igrid_query,
147  int igrid_data,
148  amrex::Real dx, amrex::Real dy, amrex::Real dz,
149  int i, int j, int k,
150  amrex::RealVect const& bcent_eb,
151  amrex::Real const state_eb,
152  amrex::Array4<amrex::Real const> const& state,
153  amrex::Array4<amrex::Real const> const& ccent,
154  amrex::Array4<amrex::EBCellFlag const> const& flag)
155 {
156  // Fitting stencil
157 
158  constexpr int dim_a = 48;
159 
160  int ii_lo = -1;
161  int jj_lo = -1;
162  int kk_lo = -1;
163  int ii_hi = 1;
164  int jj_hi = 1;
165  int kk_hi = 1;
166 
167  if (igrid_query == Vars::xvel) {
168  ii_lo = -2;
169  } else if (igrid_query == Vars::yvel) {
170  jj_lo = -2;
171  } else if (igrid_query == Vars::zvel) {
172  kk_lo = -2;
173  }
174 
175  if (igrid_data == Vars::xvel) {
176  ii_hi = 2;
177  } else if (igrid_data == Vars::yvel) {
178  jj_hi = 2;
179  } else if (igrid_data == Vars::zvel) {
180  kk_hi = 2;
181  }
182 
183  // Set bias in the index space between the two staggered grids
184  amrex::RealVect bias{0.,0.,0.};
185  bias[igrid_query-1] = 0.5;
186  bias[igrid_data-1] = -0.5;
187 
188  amrex::Real A[dim_a][AMREX_SPACEDIM];
189 
190  // Array of distances from the query point
191  int ll=0;
192  for(int kk(kk_lo); kk<=kk_hi; kk++) {
193  for(int jj(jj_lo); jj<=jj_hi; jj++) {
194  for(int ii(ii_lo); ii<=ii_hi; ii++) {
195 
196  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
197  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
198  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
199  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
200  } else {
201  A[ll][0] = amrex::Real(0.0);
202  A[ll][1] = amrex::Real(0.0);
203  A[ll][2] = amrex::Real(0.0);
204  }
205  ll++;
206  }}}
207 
208  //
209  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
210  //
211  amrex::Real du[dim_a];
212 
213  ll=0;
214  for(int kk(kk_lo); kk<=kk_hi; kk++) {
215  for(int jj(jj_lo); jj<=jj_hi; jj++) {
216  for(int ii(ii_lo); ii<=ii_hi; ii++) {
217 
218  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
219  du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
220  } else {
221  du[ll] = amrex::Real(0.0);
222  }
223  ll++;
224  }}}
225 
226  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
227  amrex::Real Atb[AMREX_SPACEDIM];
228 
229  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
230  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
231  AtA[jj][ii] = amrex::Real(0.0);
232  }
233  Atb[jj] = amrex::Real(0.0);
234  }
235 
236  for(int lc(0); lc < dim_a; ++lc)
237  {
238  AtA[0][0] += A[lc][0]* A[lc][0];
239  AtA[0][1] += A[lc][0]* A[lc][1];
240  AtA[0][2] += A[lc][0]* A[lc][2];
241  AtA[1][1] += A[lc][1]* A[lc][1];
242  AtA[1][2] += A[lc][1]* A[lc][2];
243  AtA[2][2] += A[lc][2]* A[lc][2];
244 
245  Atb[0] += A[lc][0]*du[lc];
246  Atb[1] += A[lc][1]*du[lc];
247  Atb[2] += A[lc][2]*du[lc];
248  }
249 
250  // Fill in symmetric
251  AtA[1][0] = AtA[0][1];
252  AtA[2][0] = AtA[0][2];
253  AtA[2][1] = AtA[1][2];
254 
255  amrex::Real detAtA =
256  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
257  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
258  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
259 
260  amrex::Real detAtA_x =
261  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
262  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
263  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
264 
265  // Slope at centroid of (i,j,k)
266  amrex::Real xslope = detAtA_x / detAtA;
267 
268  amrex::Real detAtA_y =
269  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
270  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
271  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
272 
273  // Slope at centroid of (i,j,k)
274  amrex::Real yslope = detAtA_y / detAtA;
275 
276  amrex::Real detAtA_z =
277  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
278  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
279  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
280 
281  // Slope at centroid of (i,j,k)
282  amrex::Real zslope = detAtA_z / detAtA;
283 
284  return {xslope,yslope,zslope};
285 }
286 
287 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
288 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
289 erf_calc_slopes_eb_staggered (int igrid_query,
290  int igrid_data,
291  amrex::Real dx, amrex::Real dy, amrex::Real dz,
292  int i, int j, int k,
293  amrex::Array4<amrex::Real const> const& state,
294  amrex::Array4<amrex::Real const> const& ccent,
295  amrex::Array4<amrex::EBCellFlag const> const& flag)
296 {
297  AMREX_ASSERT((igrid_query == Vars::xvel || igrid_query == Vars::yvel || igrid_query == Vars::zvel) && igrid_data == Vars::cons);
298 
299  // Fitting stencil
300 
301  constexpr int dim_a = 36;
302 
303  int ii_lo = -1;
304  int jj_lo = -1;
305  int kk_lo = -1;
306  int ii_hi = 1;
307  int jj_hi = 1;
308  int kk_hi = 1;
309 
310  if (igrid_query == Vars::xvel) {
311  ii_lo = -2;
312  } else if (igrid_query == Vars::yvel) {
313  jj_lo = -2;
314  } else if (igrid_query == Vars::zvel) {
315  kk_lo = -2;
316  }
317 
318  if (igrid_data == Vars::xvel) {
319  ii_hi = 2;
320  } else if (igrid_data == Vars::yvel) {
321  jj_hi = 2;
322  } else if (igrid_data == Vars::zvel) {
323  kk_hi = 2;
324  }
325 
326  // Set bias in the index space between the two staggered grids
327  amrex::RealVect bias{0.,0.,0.};
328  bias[igrid_query-1] = 0.5;
329 
330  amrex::Real A[dim_a][AMREX_SPACEDIM];
331 
332  // Array of distances from the query point
333  int ll=0;
334  for(int kk(kk_lo); kk<=kk_hi; kk++) {
335  for(int jj(jj_lo); jj<=jj_hi; jj++) {
336  for(int ii(ii_lo); ii<=ii_hi; ii++) {
337 
338  // if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
339  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
340  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - ccent(i,j,k,0) ) * dx;
341  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - ccent(i,j,k,1) ) * dy;
342  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - ccent(i,j,k,2) ) * dz;
343  } else {
344  A[ll][0] = amrex::Real(0.0);
345  A[ll][1] = amrex::Real(0.0);
346  A[ll][2] = amrex::Real(0.0);
347  }
348  ll++;
349  }}}
350 
351  //
352  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
353  //
354  amrex::Real du[dim_a];
355 
356  ll=0;
357  for(int kk(kk_lo); kk<=kk_hi; kk++) {
358  for(int jj(jj_lo); jj<=jj_hi; jj++) {
359  for(int ii(ii_lo); ii<=ii_hi; ii++) {
360 
361  // if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
362  if (!flag(i+ii,j+jj,k+kk).isCovered()) {
363  du[ll] = state(i+ii,j+jj,k+kk) - state(i,j,k);
364  } else {
365  du[ll] = amrex::Real(0.0);
366  }
367  ll++;
368  }}}
369 
370  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
371  amrex::Real Atb[AMREX_SPACEDIM];
372 
373  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
374  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
375  AtA[jj][ii] = amrex::Real(0.0);
376  }
377  Atb[jj] = amrex::Real(0.0);
378  }
379 
380  for(int lc(0); lc < dim_a; ++lc)
381  {
382  AtA[0][0] += A[lc][0]* A[lc][0];
383  AtA[0][1] += A[lc][0]* A[lc][1];
384  AtA[0][2] += A[lc][0]* A[lc][2];
385  AtA[1][1] += A[lc][1]* A[lc][1];
386  AtA[1][2] += A[lc][1]* A[lc][2];
387  AtA[2][2] += A[lc][2]* A[lc][2];
388 
389  Atb[0] += A[lc][0]*du[lc];
390  Atb[1] += A[lc][1]*du[lc];
391  Atb[2] += A[lc][2]*du[lc];
392  }
393 
394  // Fill in symmetric
395  AtA[1][0] = AtA[0][1];
396  AtA[2][0] = AtA[0][2];
397  AtA[2][1] = AtA[1][2];
398 
399  amrex::Real detAtA =
400  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
401  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
402  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
403 
404  amrex::Real detAtA_x =
405  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
406  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
407  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
408 
409  // Slope at centroid of (i,j,k)
410  amrex::Real xslope = detAtA_x / detAtA;
411 
412  amrex::Real detAtA_y =
413  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
414  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
415  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
416 
417  // Slope at centroid of (i,j,k)
418  amrex::Real yslope = detAtA_y / detAtA;
419 
420  amrex::Real detAtA_z =
421  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
422  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
423  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
424 
425  // Slope at centroid of (i,j,k)
426  amrex::Real zslope = detAtA_z / detAtA;
427 
428  return {xslope,yslope,zslope};
429 }
430 
431 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
432 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
433 erf_calc_slopes_eb_staggered_upwind (int igrid_query,
434  [[maybe_unused]] int igrid_data,
435  amrex::Real dx, amrex::Real dy, amrex::Real dz,
436  int i, int j, int k,
437  const amrex::Array4<const amrex::Real>& vel_arr,
438  amrex::Array4<amrex::Real const> const& state,
439  amrex::Array4<amrex::Real const> const& ccent,
440  amrex::Array4<amrex::EBCellFlag const> const& flag)
441 {
442  AMREX_ASSERT((igrid_query == Vars::xvel || igrid_query == Vars::yvel || igrid_query == Vars::zvel) && igrid_data == Vars::cons);
443 
444  // Fitting stencil
445 
446  int ii_lo = -1;
447  int jj_lo = -1;
448  int kk_lo = -1;
449  int ii_hi = 1;
450  int jj_hi = 1;
451  int kk_hi = 1;
452 
453  if (igrid_query == Vars::xvel && vel_arr(i,j,k)>0.0) {
454  ii_lo = -2;
455  ii_hi = 0;
456  }
457  if (igrid_query == Vars::yvel && vel_arr(i,j,k)>0.0) {
458  jj_lo = -2;
459  jj_hi = 0;
460  }
461  if (igrid_query == Vars::zvel && vel_arr(i,j,k)>0.0) {
462  kk_lo = -2;
463  kk_hi = 0;
464  }
465 
466  // Set bias in the index space between the two staggered grids
467  amrex::RealVect bias{0.,0.,0.};
468  bias[igrid_query-1] = 0.5;
469 
470  constexpr int dim_a = 27;
471 
472  amrex::Real A[dim_a][AMREX_SPACEDIM];
473 
474  // Array of distances from the query point
475  int ll=0;
476  for(int kk(kk_lo); kk<=kk_hi; kk++) {
477  for(int jj(jj_lo); jj<=jj_hi; jj++) {
478  for(int ii(ii_lo); ii<=ii_hi; ii++) {
479 
480  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
481  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - ccent(i,j,k,0) ) * dx;
482  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - ccent(i,j,k,1) ) * dy;
483  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - ccent(i,j,k,2) ) * dz;
484  } else {
485  A[ll][0] = amrex::Real(0.0);
486  A[ll][1] = amrex::Real(0.0);
487  A[ll][2] = amrex::Real(0.0);
488  }
489  ll++;
490  }}}
491 
492  //
493  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
494  //
495  amrex::Real du[dim_a];
496 
497  ll=0;
498  for(int kk(kk_lo); kk<=kk_hi; kk++) {
499  for(int jj(jj_lo); jj<=jj_hi; jj++) {
500  for(int ii(ii_lo); ii<=ii_hi; ii++) {
501 
502  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
503  du[ll] = state(i+ii,j+jj,k+kk) - state(i,j,k);
504  } else {
505  du[ll] = amrex::Real(0.0);
506  }
507  ll++;
508  }}}
509 
510  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
511  amrex::Real Atb[AMREX_SPACEDIM];
512 
513  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
514  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
515  AtA[jj][ii] = amrex::Real(0.0);
516  }
517  Atb[jj] = amrex::Real(0.0);
518  }
519 
520  for(int lc(0); lc < dim_a; ++lc)
521  {
522  AtA[0][0] += A[lc][0]* A[lc][0];
523  AtA[0][1] += A[lc][0]* A[lc][1];
524  AtA[0][2] += A[lc][0]* A[lc][2];
525  AtA[1][1] += A[lc][1]* A[lc][1];
526  AtA[1][2] += A[lc][1]* A[lc][2];
527  AtA[2][2] += A[lc][2]* A[lc][2];
528 
529  Atb[0] += A[lc][0]*du[lc];
530  Atb[1] += A[lc][1]*du[lc];
531  Atb[2] += A[lc][2]*du[lc];
532  }
533 
534  // Fill in symmetric
535  AtA[1][0] = AtA[0][1];
536  AtA[2][0] = AtA[0][2];
537  AtA[2][1] = AtA[1][2];
538 
539  amrex::Real detAtA =
540  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
541  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
542  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
543 
544  amrex::Real detAtA_x =
545  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
546  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
547  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
548 
549  // Slope at centroid of (i,j,k)
550  amrex::Real xslope = detAtA_x / detAtA;
551 
552  amrex::Real detAtA_y =
553  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
554  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
555  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
556 
557  // Slope at centroid of (i,j,k)
558  amrex::Real yslope = detAtA_y / detAtA;
559 
560  amrex::Real detAtA_z =
561  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
562  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
563  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
564 
565  // Slope at centroid of (i,j,k)
566  amrex::Real zslope = detAtA_z / detAtA;
567 
568  return {xslope,yslope,zslope};
569 }
570 
571 #endif
erf_calc_slopes_eb_Dirichlet
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet(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)
Definition: ERF_EBSlopes.H:20
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)
Definition: ERF_EBSlopes.H:146
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)
Definition: ERF_EBSlopes.H:433
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)
Definition: ERF_EBSlopes.H:289
ERF_IndexDefines.H
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::cons
@ cons
Definition: ERF_IndexDefines.H:140
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:143
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142