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