Block-Structured AMR Software Framework
AMReX_Habec_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_HABEC_3D_H_
2 #define AMREX_HABEC_3D_H_
3 #include <AMReX_Config.H>
4 
5 #ifdef AMREX_USE_EB
6 #include <AMReX_EBMultiFabUtil.H>
7 #include <AMReX_MultiCutFab.H>
8 #include <AMReX_EBFabFactory.H>
9 #include <AMReX_MLLinOp_K.H>
10 #endif
11 
12 namespace amrex {
13 
14 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
15 void habec_mat (GpuArray<Real,2*AMREX_SPACEDIM+1>& sten, int i, int j, int k,
16  Dim3 const& boxlo, Dim3 const& boxhi,
17  Real sa, Array4<Real const> const& a,
18  Real sb, GpuArray<Real,AMREX_SPACEDIM> const& dx,
19  GpuArray<Array4<Real const>, AMREX_SPACEDIM> const& b,
20  GpuArray<int,AMREX_SPACEDIM*2> const& bctype,
21  GpuArray<Real,AMREX_SPACEDIM*2> const& bcl, int bho,
22  GpuArray<Array4<int const>, AMREX_SPACEDIM*2> const& msk,
23  Array4<Real> const& diaginv)
24 {
25  sten[1] = -(sb / (dx[0]*dx[0])) * b[0](i,j,k);
26  sten[2] = -(sb / (dx[0]*dx[0])) * b[0](i+1,j,k);
27  sten[3] = -(sb / (dx[1]*dx[1])) * b[1](i,j,k);
28  sten[4] = -(sb / (dx[1]*dx[1])) * b[1](i,j+1,k);
29  sten[5] = -(sb / (dx[2]*dx[2])) * b[2](i,j,k);
30  sten[6] = -(sb / (dx[2]*dx[2])) * b[2](i,j,k+1);
31  sten[0] = -(sten[1] + sten[2] + sten[3] + sten[4] + sten[5] + sten[6]);
32  if (sa != Real(0.0)) {
33  sten[0] += sa * a(i,j,k);
34  }
35 
36  // xlo
37  if (i == boxlo.x) {
38  int cdir = Orientation(Direction::x, Orientation::low);
39  if (msk[cdir](i-1,j,k) > 0) {
40  Real bf1, bf2;
41  detail::comp_bf(bf1, bf2, sb, dx[0], bctype[cdir], bcl[cdir], bho);
42  sten[0] += bf1 * b[0](i,j,k);
43  sten[1] = Real(0.0);
44  sten[2] += bf2 * b[0](i,j,k);
45  }
46  }
47 
48  // xhi
49  if (i == boxhi.x) {
50  int cdir = Orientation(Direction::x, Orientation::high);
51  if (msk[cdir](i+1,j,k) > 0) {
52  Real bf1, bf2;
53  detail::comp_bf(bf1, bf2, sb, dx[0], bctype[cdir], bcl[cdir], bho);
54  sten[0] += bf1 * b[0](i+1,j,k);
55  sten[1] += bf2 * b[0](i+1,j,k);
56  sten[2] = Real(0.0);
57  }
58  }
59 
60  // ylo
61  if (j == boxlo.y) {
62  int cdir = Orientation(Direction::y, Orientation::low);
63  if (msk[cdir](i,j-1,k) > 0) {
64  Real bf1, bf2;
65  detail::comp_bf(bf1, bf2, sb, dx[1], bctype[cdir], bcl[cdir], bho);
66  sten[0] += bf1 * b[1](i,j,k);
67  sten[3] = Real(0.0);
68  sten[4] += bf2 * b[1](i,j,k);
69  }
70  }
71 
72  // yhi
73  if (j == boxhi.y) {
74  int cdir = Orientation(Direction::y, Orientation::high);
75  if (msk[cdir](i,j+1,k) > 0) {
76  Real bf1, bf2;
77  detail::comp_bf(bf1, bf2, sb, dx[1], bctype[cdir], bcl[cdir], bho);
78  sten[0] += bf1 * b[1](i,j+1,k);
79  sten[3] += bf2 * b[1](i,j+1,k);
80  sten[4] = Real(0.0);
81  }
82  }
83 
84  // zlo
85  if (k == boxlo.z) {
86  int cdir = Orientation(Direction::z, Orientation::low);
87  if (msk[cdir](i,j,k-1) > 0) {
88  Real bf1, bf2;
89  detail::comp_bf(bf1, bf2, sb, dx[2], bctype[cdir], bcl[cdir], bho);
90  sten[0] += bf1 * b[2](i,j,k);
91  sten[5] = Real(0.0);
92  sten[6] += bf2 * b[2](i,j,k);
93  }
94  }
95 
96  // zhi
97  if (k == boxhi.z) {
98  int cdir = Orientation(Direction::z, Orientation::high);
99  if (msk[cdir](i,j,k+1) > 0) {
100  Real bf1, bf2;
101  detail::comp_bf(bf1, bf2, sb, dx[2], bctype[cdir], bcl[cdir], bho);
102  sten[0] += bf1 * b[2](i,j,k+1);
103  sten[5] += bf2 * b[2](i,j,k+1);
104  sten[6] = Real(0.0);
105  }
106  }
107 
108  diaginv(i,j,k) = Real(1.0) / sten[0];
109  sten[0] = Real(1.0);
110  for (int m = 1; m < 2*AMREX_SPACEDIM+1; ++m) {
111  sten[m] *= diaginv(i,j,k);
112  }
113 }
114 
115 template <typename Int>
116 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
117 void habec_ijmat (GpuArray<Real,2*AMREX_SPACEDIM+1>& sten, Array4<Int> const& ncols,
118  Array4<Real> const& diaginv, int i, int j, int k,
119  Array4<Int const> const& cell_id,
120  Real sa, Array4<Real const> const& a,
121  Real sb, GpuArray<Real,AMREX_SPACEDIM> const& dx,
122  GpuArray<Array4<Real const>, AMREX_SPACEDIM> const& b,
123  GpuArray<int,AMREX_SPACEDIM*2> const& bctype,
124  GpuArray<Real,AMREX_SPACEDIM*2> const& bcl, int bho,
125  Array4<int const> const& osm)
126 {
127  if (!osm || osm(i,j,k) != 0) {
128  sten[1] = -(sb / (dx[0]*dx[0])) * b[0](i,j,k);
129  sten[2] = -(sb / (dx[0]*dx[0])) * b[0](i+1,j,k);
130  sten[3] = -(sb / (dx[1]*dx[1])) * b[1](i,j,k);
131  sten[4] = -(sb / (dx[1]*dx[1])) * b[1](i,j+1,k);
132  sten[5] = -(sb / (dx[2]*dx[2])) * b[2](i,j,k);
133  sten[6] = -(sb / (dx[2]*dx[2])) * b[2](i,j,k+1);
134  sten[0] = -(sten[1] + sten[2] + sten[3] + sten[4] + sten[5] + sten[6]);
135  if (sa != Real(0.0)) {
136  sten[0] += sa * a(i,j,k);
137  }
138 
139  // xlo
140  if (cell_id(i-1,j,k) < 0) {
141  int cdir = Orientation(Direction::x, Orientation::low);
142  Real bf1, bf2;
143  detail::comp_bf(bf1, bf2, sb, dx[0], bctype[cdir], bcl[cdir], bho);
144  sten[0] += bf1 * b[0](i,j,k);
145  sten[1] = Real(0.0);
146  sten[2] += bf2 * b[0](i,j,k);
147  }
148 
149  // xhi
150  if (cell_id(i+1,j,k) < 0) {
151  int cdir = Orientation(Direction::x, Orientation::high);
152  Real bf1, bf2;
153  detail::comp_bf(bf1, bf2, sb, dx[0], bctype[cdir], bcl[cdir], bho);
154  sten[0] += bf1 * b[0](i+1,j,k);
155  sten[1] += bf2 * b[0](i+1,j,k);
156  sten[2] = Real(0.0);
157  }
158 
159  // ylo
160  if (cell_id(i,j-1,k) < 0) {
161  int cdir = Orientation(Direction::y, Orientation::low);
162  Real bf1, bf2;
163  detail::comp_bf(bf1, bf2, sb, dx[1], bctype[cdir], bcl[cdir], bho);
164  sten[0] += bf1 * b[1](i,j,k);
165  sten[3] = Real(0.0);
166  sten[4] += bf2 * b[1](i,j,k);
167  }
168 
169  // yhi
170  if (cell_id(i,j+1,k) < 0) {
171  int cdir = Orientation(Direction::y, Orientation::high);
172  Real bf1, bf2;
173  detail::comp_bf(bf1, bf2, sb, dx[1], bctype[cdir], bcl[cdir], bho);
174  sten[0] += bf1 * b[1](i,j+1,k);
175  sten[3] += bf2 * b[1](i,j+1,k);
176  sten[4] = Real(0.0);
177  }
178 
179  // zlo
180  if (cell_id(i,j,k-1) < 0) {
181  int cdir = Orientation(Direction::z, Orientation::low);
182  Real bf1, bf2;
183  detail::comp_bf(bf1, bf2, sb, dx[2], bctype[cdir], bcl[cdir], bho);
184  sten[0] += bf1 * b[2](i,j,k);
185  sten[5] = Real(0.0);
186  sten[6] += bf2 * b[2](i,j,k);
187  }
188 
189  // zhi
190  if (cell_id(i,j,k+1) < 0) {
191  int cdir = Orientation(Direction::z, Orientation::high);
192  Real bf1, bf2;
193  detail::comp_bf(bf1, bf2, sb, dx[2], bctype[cdir], bcl[cdir], bho);
194  sten[0] += bf1 * b[2](i,j,k+1);
195  sten[5] += bf2 * b[2](i,j,k+1);
196  sten[6] = Real(0.0);
197  }
198  } else {
199  sten[0] = Real(1.0);
200  for (int m = 1; m < 2*AMREX_SPACEDIM+1; ++m) {
201  sten[m] = Real(0.0);
202  }
203  }
204 
205  diaginv(i,j,k) = Real(1.0) / sten[0];
206  sten[0] = Real(1.0);
207  for (int m = 1; m < 2*AMREX_SPACEDIM+1; ++m) {
208  sten[m] *= diaginv(i,j,k);
209  }
210 
211  ncols(i,j,k) = 0;
212  for (int m = 0; m < 2*AMREX_SPACEDIM+1; ++m) {
213  ncols(i,j,k) += (sten[m] != Real(0.0));
214  }
215 }
216 
217 template <typename Int>
218 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
219 void habec_cols (GpuArray<Int,2*AMREX_SPACEDIM+1>& sten, int i, int j, int k,
220  Array4<Int const> const& cell_id)
221 {
222  sten[0] = cell_id(i ,j ,k );
223  sten[1] = cell_id(i-1,j ,k );
224  sten[2] = cell_id(i+1,j ,k );
225  sten[3] = cell_id(i ,j-1,k );
226  sten[4] = cell_id(i ,j+1,k );
227  sten[5] = cell_id(i ,j ,k-1);
228  sten[6] = cell_id(i ,j ,k+1);
229 }
230 
231 #ifdef AMREX_USE_EB
232 
233 template <typename Int>
234 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
235 void habec_cols_eb (GpuArray<Int,27>& sten, int i, int j, int k,
236  Array4<Int const> const& cell_id)
237 {
238  sten[0 ] = cell_id(i-1,j-1,k-1);
239  sten[1 ] = cell_id(i ,j-1,k-1);
240  sten[2 ] = cell_id(i+1,j-1,k-1);
241  sten[3 ] = cell_id(i-1,j ,k-1);
242  sten[4 ] = cell_id(i ,j ,k-1);
243  sten[5 ] = cell_id(i+1,j ,k-1);
244  sten[6 ] = cell_id(i-1,j+1,k-1);
245  sten[7 ] = cell_id(i ,j+1,k-1);
246  sten[8 ] = cell_id(i+1,j+1,k-1);
247  sten[9 ] = cell_id(i-1,j-1,k );
248  sten[10] = cell_id(i ,j-1,k );
249  sten[11] = cell_id(i+1,j-1,k );
250  sten[12] = cell_id(i-1,j ,k );
251  sten[13] = cell_id(i ,j ,k );
252  sten[14] = cell_id(i+1,j ,k );
253  sten[15] = cell_id(i-1,j+1,k );
254  sten[16] = cell_id(i ,j+1,k );
255  sten[17] = cell_id(i+1,j+1,k );
256  sten[18] = cell_id(i-1,j-1,k+1);
257  sten[19] = cell_id(i ,j-1,k+1);
258  sten[20] = cell_id(i+1,j-1,k+1);
259  sten[21] = cell_id(i-1,j ,k+1);
260  sten[22] = cell_id(i ,j ,k+1);
261  sten[23] = cell_id(i+1,j ,k+1);
262  sten[24] = cell_id(i-1,j+1,k+1);
263  sten[25] = cell_id(i ,j+1,k+1);
264  sten[26] = cell_id(i+1,j+1,k+1);
265 }
266 
267 template <typename Int>
268 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
269 void habec_ijmat_eb (GpuArray<Real,27>& sten, Array4<Int> const& ncols,
270  Array4<Real> const& diaginv, int i, int j, int k,
271  Array4<Int const> const& cell_id,
272  Real sa, Array4<Real const> const& a,
273  Real sb, GpuArray<Real,AMREX_SPACEDIM> const& dx,
274  GpuArray<Array4<Real const>, AMREX_SPACEDIM> const& b,
275  GpuArray<int,AMREX_SPACEDIM*2> const& bctype,
276  GpuArray<Real,AMREX_SPACEDIM*2> const& bcl, int bho,
277  Array4<const EBCellFlag> const& flag,
278  Array4<Real const> const& vfrc,
279  Array4<Real const> const& apx,
280  Array4<Real const> const& apy,
281  Array4<Real const> const& apz,
282  Array4<Real const> const& fcx,
283  Array4<Real const> const& fcy,
284  Array4<Real const> const& fcz,
285  Array4<Real const> const& ba,
286  Array4<Real const> const& bcen,
287  Array4<Real const> const& beb)
288 {
289  for (int m = 0; m < 27; ++m) {
290  sten[m] = Real(0.0);
291  }
292 
293  auto& mat_tmp = reinterpret_cast<Array3D<Real,-1,1,-1,1,-1,1>&>(sten);
294 
295  GpuArray<Real,AMREX_SPACEDIM> fac{sb/(dx[0]*dx[0]), sb/(dx[1]*dx[1]), sb/(dx[2]*dx[2])};
296 
297  if (flag(i,j,k).isRegular())
298  {
299  mat_tmp(0,0,0) = fac[0]*(b[0](i,j,k)+b[0](i+1,j,k))
300  + fac[1]*(b[1](i,j,k)+b[1](i,j+1,k))
301  + fac[2]*(b[2](i,j,k)+b[2](i,j,k+1));
302  mat_tmp(-1, 0, 0) = -fac[0]*b[0](i,j,k);
303  mat_tmp( 1, 0, 0) = -fac[0]*b[0](i+1,j,k);
304  mat_tmp( 0,-1, 0) = -fac[1]*b[1](i,j,k);
305  mat_tmp( 0, 1, 0) = -fac[1]*b[1](i,j+1,k);
306  mat_tmp( 0, 0,-1) = -fac[2]*b[2](i,j,k);
307  mat_tmp( 0, 0, 1) = -fac[2]*b[2](i,j,k+1);
308 
309  if (cell_id(i-1,j,k) < 0) {
310  int cdir = Orientation(Direction::x, Orientation::low);
311  Real bf1, bf2;
312  detail::comp_bf(bf1, bf2, sb, dx[0], bctype[cdir], bcl[cdir], bho);
313  mat_tmp(0,0,0) += bf1*b[0](i,j,k);
314  mat_tmp(-1,0,0) = Real(0.0);
315  mat_tmp(1,0,0) += bf2*b[0](i,j,k);
316  }
317 
318  if (cell_id(i+1,j,k) < 0) {
319  int cdir = Orientation(Direction::x, Orientation::high);
320  Real bf1, bf2;
321  detail::comp_bf(bf1, bf2, sb, dx[0], bctype[cdir], bcl[cdir], bho);
322  mat_tmp(0,0,0) += bf1*b[0](i+1,j,k);
323  mat_tmp(1,0,0) = Real(0.0);
324  mat_tmp(-1,0,0) += bf2*b[0](i+1,j,k);
325  }
326 
327  if (cell_id(i,j-1,k) < 0) {
328  int cdir = Orientation(Direction::y, Orientation::low);
329  Real bf1, bf2;
330  detail::comp_bf(bf1, bf2, sb, dx[1], bctype[cdir], bcl[cdir], bho);
331  mat_tmp(0,0,0) += bf1*b[1](i,j,k);
332  mat_tmp(0,-1,0) = Real(0.0);
333  mat_tmp(0,1,0) += bf2*b[1](i,j,k);
334  }
335 
336  if (cell_id(i,j+1,k) < 0) {
337  int cdir = Orientation(Direction::y, Orientation::high);
338  Real bf1, bf2;
339  detail::comp_bf(bf1, bf2, sb, dx[1], bctype[cdir], bcl[cdir], bho);
340  mat_tmp(0,0,0) += bf1*b[1](i,j+1,k);
341  mat_tmp(0,1,0) = Real(0.0);
342  mat_tmp(0,-1,0) += bf2*b[1](i,j+1,k);
343  }
344 
345  if (cell_id(i,j,k-1) < 0) {
346  int cdir = Orientation(Direction::z, Orientation::low);
347  Real bf1, bf2;
348  detail::comp_bf(bf1, bf2, sb, dx[2], bctype[cdir], bcl[cdir], bho);
349  mat_tmp(0,0,0) += bf1*b[2](i,j,k);
350  mat_tmp(0,0,-1) = Real(0.0);
351  mat_tmp(0,0,1) += bf2*b[2](i,j,k);
352  }
353 
354  if (cell_id(i,j,k+1) < 0) {
355  int cdir = Orientation(Direction::z, Orientation::high);
356  Real bf1, bf2;
357  detail::comp_bf(bf1, bf2, sb, dx[2], bctype[cdir], bcl[cdir], bho);
358  mat_tmp(0,0,0) += bf1*b[2](i,j,k+1);
359  mat_tmp(0,0,1) = Real(0.0);
360  mat_tmp(0,0,-1) += bf2*b[2](i,j,k+1);
361  }
362 
363  if (sa != Real(0.0)) {
364  mat_tmp(0,0,0) += sa*a(i,j,k);
365  }
366  }
367  else if (flag(i,j,k).isSingleValued())
368  {
369  int cdir = Orientation(Direction::x, Orientation::low);
370  Real area = apx(i,j,k);
371 
372  if (area > Real(0.0))
373  {
374  int joff, koff, jj, kk;
375  Real fracy, fracz;
376  Real f = fac[0];
377 
378  if (area != Real(1.0)) {
379  joff = static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,0)));
380  koff = static_cast<int>(std::copysign(Real(1.0), fcx(i,j,k,1)));
381  jj = j+joff;
382  kk = k+koff;
383  if (cell_id(i-1,jj,k) < 0 && cell_id(i,jj,k) < 0) {
384  fracy = Real(0.0);
385  } else {
386  fracy = std::abs(fcx(i,j,k,0));
387  }
388  if (cell_id(i-1,j,kk) < 0 && cell_id(i,j,kk) < 0) {
389  fracz = Real(0.0);
390  } else {
391  fracz = std::abs(fcx(i,j,k,1));
392  }
393  if (cell_id(i-1,jj,kk) < 0 && cell_id(i,jj,kk) < 0 && (fracy*fracz) > Real(0.0)) {
394  fracy = Real(0.0);
395  fracz = Real(0.0);
396  }
397  } else {
398  joff = 0;
399  koff = 0;
400  jj = j;
401  kk = k;
402  fracy = Real(0.0);
403  fracz = Real(0.0);
404  }
405 
406  Real bf1, bf2, bflo;
407  detail::comp_bflo(bf1, bf2, bflo, sb, dx[0], bctype[cdir], bcl[cdir], bho);
408 
409  Real tmp = (Real(1.0)-fracy)*(1.0-fracz)*area*b[0](i,j,k);
410  // cell(i-1,j,k) is not covered because area > 0
411  if (cell_id(i-1,j,k) >= 0) {
412  mat_tmp( 0,0,0) += tmp*f;
413  mat_tmp(-1,0,0) -= tmp*f;
414  } else if (cell_id(i+1,j,k) < 0 || apx(i+1,j,k) == Real(0.0)) {
415  mat_tmp(0,0,0) += tmp*(f+bflo);
416  } else {
417  mat_tmp(0,0,0) += tmp*(f+bf1);
418  mat_tmp(1,0,0) += tmp* bf2;
419  }
420 
421  if (fracy > Real(0.0)) {
422  tmp = fracy*(Real(1.0)-fracz)*area*b[0](i,jj,k);
423  if (cell_id(i-1,jj,k) >= 0 && cell_id(i,jj,k) >= 0) {
424  mat_tmp(-1,joff,0) -= tmp*f;
425  mat_tmp( 0,joff,0) += tmp*f;
426  } else if (cell_id(i+1,jj,k) < 0 || apx(i+1,jj,k) == Real(0.0)) {
427  mat_tmp( 0,joff,0) += tmp*(f+bflo);
428  } else {
429  mat_tmp( 0,joff,0) += tmp*(f+bf1);
430  mat_tmp( 1,joff,0) += tmp* bf2;
431  }
432  }
433 
434  if (fracz > Real(0.0)) {
435  tmp = fracz*(Real(1.0)-fracy)*area*b[0](i,j,kk);
436  if (cell_id(i-1,j,kk) >= 0 && cell_id(i,j,kk) >= 0) {
437  mat_tmp(-1,0,koff) -= tmp*f;
438  mat_tmp( 0,0,koff) += tmp*f;
439  } else if (cell_id(i+1,j,kk) < 0 || apx(i+1,j,kk) == Real(0.0)) {
440  mat_tmp( 0,0,koff) += tmp*(f+bflo);
441  } else {
442  mat_tmp( 0,0,koff) += tmp*(f+bf1);
443  mat_tmp( 1,0,koff) += tmp* bf2;
444  }
445  }
446 
447  if (fracy > Real(0.0) && fracz > Real(0.0)) {
448  tmp = fracy*fracz*area*b[0](i,jj,kk);
449  if (cell_id(i-1,jj,kk) >= 0 && cell_id(i,jj,kk) >= 0) {
450  mat_tmp(-1,joff,koff) -= tmp*f;
451  mat_tmp( 0,joff,koff) += tmp*f;
452  } else if (cell_id(i+1,jj,kk) < 0 || apx(i+1,jj,kk) == Real(0.0)) {
453  mat_tmp( 0,joff,koff) += tmp*(f+bflo);
454  } else {
455  mat_tmp( 0,joff,koff) += tmp*(f+bf1);
456  mat_tmp( 1,joff,koff) += tmp* bf2;
457  }
458  }
459  }
460 
461  cdir = Orientation(Direction::x, Orientation::high);
462  area = apx(i+1,j,k);
463 
464  if (area > Real(0.0))
465  {
466  int joff, koff, jj, kk;
467  Real fracy, fracz;
468  Real f = fac[0];
469 
470  if (area != Real(1.0)) {
471  joff = static_cast<int>(std::copysign(Real(1.0), fcx(i+1,j,k,0)));
472  koff = static_cast<int>(std::copysign(Real(1.0), fcx(i+1,j,k,1)));
473  jj = j+joff;
474  kk = k+koff;
475  if (cell_id(i,jj,k) < 0 && cell_id(i+1,jj,k) < 0) {
476  fracy = Real(0.0);
477  } else {
478  fracy = std::abs(fcx(i+1,j,k,0));
479  }
480  if (cell_id(i,j,kk) < 0 && cell_id(i+1,j,kk) < 0) {
481  fracz = Real(0.0);
482  } else {
483  fracz = std::abs(fcx(i+1,j,k,1));
484  }
485  if (cell_id(i,jj,kk) < 0 && cell_id(i+1,jj,kk) < 0 && (fracy*fracz) > Real(0.0)) {
486  fracy = Real(0.0);
487  fracz = Real(0.0);
488  }
489  } else {
490  joff = 0;
491  koff = 0;
492  jj = j;
493  kk = k;
494  fracy = Real(0.0);
495  fracz = Real(0.0);
496  }
497 
498  Real bf1, bf2, bflo;
499  detail::comp_bflo(bf1, bf2, bflo, sb, dx[0], bctype[cdir], bcl[cdir], bho);
500 
501  Real tmp = (Real(1.0)-fracy)*(Real(1.0)-fracz)*area*b[0](i+1,j,k);
502  if (cell_id(i+1,j,k) >= 0) {
503  mat_tmp(0,0,0) += tmp*f;
504  mat_tmp(1,0,0) -= tmp*f;
505  } else if (cell_id(i-1,j,k) < 0 || apx(i-1,j,k) == Real(0.0)) {
506  mat_tmp(0,0,0) += tmp*(f+bflo);
507  } else {
508  mat_tmp( 0,0,0) += tmp*(f+bf1);
509  mat_tmp(-1,0,0) += tmp* bf2;
510  }
511 
512  if (fracy > Real(0.0)) {
513  tmp = fracy*(Real(1.0)-fracz)*area*b[0](i+1,jj,k);
514  if (cell_id(i,jj,k) >= 0 && cell_id(i+1,jj,k) >= 0) {
515  mat_tmp(0,joff,0) += tmp*f;
516  mat_tmp(1,joff,0) -= tmp*f;
517  } else if (cell_id(i-1,jj,k) < 0 || apx(i-1,jj,k) == Real(0.0)) {
518  mat_tmp(0,joff,0) += tmp*(f+bflo);
519  } else {
520  mat_tmp( 0,joff,0) += tmp*(f+bf1);
521  mat_tmp(-1,joff,0) += tmp* bf2;
522  }
523  }
524 
525  if (fracz > Real(0.0)) {
526  tmp = fracz*(Real(1.0)-fracy)*area*b[0](i+1,j,kk);
527  if (cell_id(i,j,kk) >= 0 && cell_id(i+1,j,kk) >= 0) {
528  mat_tmp(0,0,koff) += tmp*f;
529  mat_tmp(1,0,koff) -= tmp*f;
530  } else if (cell_id(i-1,j,kk) < 0 || apx(i-1,j,kk) == Real(0.0)) {
531  mat_tmp(0,0,koff) += tmp*(f+bflo);
532  } else {
533  mat_tmp( 0,0,koff) += tmp*(f+bf1);
534  mat_tmp(-1,0,koff) += tmp* bf2;
535  }
536  }
537 
538  if (fracy > Real(0.0) && fracz > Real(0.0)) {
539  tmp = fracy*fracz*area*b[0](i+1,jj,kk);
540  if (cell_id(i,jj,kk) >= 0 && cell_id(i+1,jj,kk) >= 0) {
541  mat_tmp(0,joff,koff) += tmp*f;
542  mat_tmp(1,joff,koff) -= tmp*f;
543  } else if (cell_id(i-1,jj,kk) < 0 || apx(i-1,jj,kk) == Real(0.0)) {
544  mat_tmp(0,joff,koff) += tmp*(f+bflo);
545  } else {
546  mat_tmp( 0,joff,koff) += tmp*(f+bf1);
547  mat_tmp(-1,joff,koff) += tmp* bf2;
548  }
549  }
550  }
551 
552  cdir = Orientation(Direction::y, Orientation::low);
553  area = apy(i,j,k);
554 
555  if (area > Real(0.0))
556  {
557  int ioff, koff, ii, kk;
558  Real fracx, fracz;
559  Real f = fac[1];
560 
561  if (area != Real(1.0)) {
562  ioff = static_cast<int>(std::copysign(Real(1.0), fcy(i,j,k,0)));
563  koff = static_cast<int>(std::copysign(Real(1.0), fcy(i,j,k,1)));
564  ii = i+ioff;
565  kk = k+koff;
566  if (cell_id(ii,j-1,k) < 0 && cell_id(ii,j,k) < 0) {
567  fracx = Real(0.0);
568  } else {
569  fracx = std::abs(fcy(i,j,k,0));
570  }
571  if (cell_id(i,j-1,kk) < 0 && cell_id(i,j,kk) < 0) {
572  fracz = Real(0.0);
573  } else {
574  fracz = std::abs(fcy(i,j,k,1));
575  }
576  if (cell_id(ii,j-1,kk) < 0 && cell_id(ii,j,kk) < 0 && fracx*fracz > Real(0.0)) {
577  fracx = Real(0.0);
578  fracz = Real(0.0);
579  }
580  } else {
581  ioff = 0;
582  koff = 0;
583  ii = i;
584  kk = k;
585  fracx = Real(0.0);
586  fracz = Real(0.0);
587  }
588 
589  Real bf1, bf2, bflo;
590  detail::comp_bflo(bf1, bf2, bflo, sb, dx[1], bctype[cdir], bcl[cdir], bho);
591 
592  Real tmp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*area*b[1](i,j,k);
593  if (cell_id(i,j-1,k) >= 0) {
594  mat_tmp(0, 0,0) += tmp*f;
595  mat_tmp(0,-1,0) -= tmp*f;
596  } else if (cell_id(i,j+1,k) < 0 || apy(i,j+1,k) == Real(0.0)) {
597  mat_tmp(0,0,0) += tmp*(f+bflo);
598  } else {
599  mat_tmp(0,0,0) += tmp*(f+bf1);
600  mat_tmp(0,1,0) += tmp* bf2;
601  }
602 
603  if (fracx > Real(0.0)) {
604  tmp = fracx*(Real(1.0)-fracz)*area*b[1](ii,j,k);
605  if (cell_id(ii,j-1,k) >= 0 && cell_id(ii,j,k) >= 0) {
606  mat_tmp(ioff,-1,0) -= tmp*f;
607  mat_tmp(ioff, 0,0) += tmp*f;
608  } else if (cell_id(ii,j+1,k) < 0 || apy(ii,j+1,k) == Real(0.0)) {
609  mat_tmp(ioff,0,0) += tmp*(f+bflo);
610  } else {
611  mat_tmp(ioff,0,0) += tmp*(f+bf1);
612  mat_tmp(ioff,1,0) += tmp* bf2;
613  }
614  }
615 
616  if (fracz > Real(0.0)) {
617  tmp = fracz*(Real(1.0)-fracx)*area*b[1](i,j,kk);
618  if (cell_id(i,j-1,kk) >= 0 && cell_id(i,j,kk) >= 0) {
619  mat_tmp(0,-1,koff) -= tmp*f;
620  mat_tmp(0, 0,koff) += tmp*f;
621  } else if (cell_id(i,j+1,kk) < 0 || apy(i,j+1,kk) == Real(0.0)) {
622  mat_tmp(0,0,koff) += tmp*(f+bflo);
623  } else {
624  mat_tmp(0,0,koff) += tmp*(f+bf1);
625  mat_tmp(0,1,koff) += tmp* bf2;
626  }
627  }
628 
629  if (fracx > Real(0.0) && fracz > Real(0.0)) {
630  tmp = fracx*fracz*area*b[1](ii,j,kk);
631  if (cell_id(ii,j-1,kk) >= 0 && cell_id(ii,j,kk) >= 0) {
632  mat_tmp(ioff,-1,koff) -= tmp*f;
633  mat_tmp(ioff, 0,koff) += tmp*f;
634  } else if (cell_id(ii,j+1,kk) < 0 || apy(ii,j+1,kk) == Real(0.0)) {
635  mat_tmp(ioff,0,koff) += tmp*(f+bflo);
636  } else {
637  mat_tmp(ioff,0,koff) += tmp*(f+bf1);
638  mat_tmp(ioff,1,koff) += tmp* bf2;
639  }
640  }
641  }
642 
643  cdir = Orientation(Direction::y, Orientation::high);
644  area = apy(i,j+1,k);
645 
646  if (area > Real(0.0))
647  {
648  int ioff, koff, ii, kk;
649  Real fracx, fracz;
650  Real f = fac[1];
651 
652  if (area != Real(1.0)) {
653  ioff = static_cast<int>(std::copysign(Real(1.0), fcy(i,j+1,k,0)));
654  koff = static_cast<int>(std::copysign(Real(1.0), fcy(i,j+1,k,1)));
655  ii = i+ioff;
656  kk = k+koff;
657  if (cell_id(ii,j,k) < 0 && cell_id(ii,j+1,k) < 0) {
658  fracx = Real(0.0);
659  } else {
660  fracx = std::abs(fcy(i,j+1,k,0));
661  }
662  if (cell_id(i,j,kk) < 0 && cell_id(i,j+1,kk) < 0) {
663  fracz = Real(0.0);
664  } else {
665  fracz = std::abs(fcy(i,j+1,k,1));
666  }
667  if (cell_id(ii,j,kk) < 0 && cell_id(ii,j+1,kk) < 0 && fracx*fracz > Real(0.0)) {
668  fracx = Real(0.0);
669  fracz = Real(0.0);
670  }
671  } else {
672  ioff = 0;
673  koff = 0;
674  ii = i;
675  kk = k;
676  fracx = Real(0.0);
677  fracz = Real(0.0);
678  }
679 
680  Real bf1, bf2, bflo;
681  detail::comp_bflo(bf1, bf2, bflo, sb, dx[1], bctype[cdir], bcl[cdir], bho);
682 
683  Real tmp = (Real(1.0)-fracx)*(Real(1.0)-fracz)*area*b[1](i,j+1,k);
684  if (cell_id(i,j+1,k) >= 0) {
685  mat_tmp(0,0,0) += tmp*f;
686  mat_tmp(0,1,0) -= tmp*f;
687  } else if (cell_id(i,j-1,k) < 0 || apy(i,j-1,k) == Real(0.0)) {
688  mat_tmp(0,0,0) += tmp*(f+bflo);
689  } else {
690  mat_tmp(0, 0,0) += tmp*(f+bf1);
691  mat_tmp(0,-1,0) += tmp* bf2;
692  }
693 
694  if (fracx > Real(0.0)) {
695  tmp = fracx*(Real(1.0)-fracz)*area*b[1](ii,j+1,k);
696  if (cell_id(ii,j,k) >= 0 && cell_id(ii,j+1,k) >= 0) {
697  mat_tmp(ioff,0,0) += tmp*f;
698  mat_tmp(ioff,1,0) -= tmp*f;
699  } else if (cell_id(ii,j-1,k) < 0 || apy(ii,j-1,k) == Real(0.0)) {
700  mat_tmp(ioff,0,0) += tmp*(f+bflo);
701  } else {
702  mat_tmp(ioff, 0,0) += tmp*(f+bf1);
703  mat_tmp(ioff,-1,0) += tmp* bf2;
704  }
705  }
706 
707  if (fracz > Real(0.0)) {
708  tmp = fracz*(Real(1.0)-fracx)*area*b[1](i,j+1,kk);
709  if (cell_id(i,j,kk) >= 0 && cell_id(i,j+1,kk) >= 0) {
710  mat_tmp(0,0,koff) += tmp*f;
711  mat_tmp(0,1,koff) -= tmp*f;
712  } else if (cell_id(i,j-1,kk) < 0 || apy(i,j-1,kk) == Real(0.0)) {
713  mat_tmp(0,0,koff) += tmp*(f+bflo);
714  } else {
715  mat_tmp(0, 0,koff) += tmp*(f+bf1);
716  mat_tmp(0,-1,koff) += tmp* bf2;
717  }
718  }
719 
720  if (fracx > Real(0.0) && fracz > Real(0.0)) {
721  tmp = fracx*fracz*area*b[1](ii,j+1,kk);
722  if (cell_id(ii,j,kk) >= 0 && cell_id(ii,j+1,kk) >= 0) {
723  mat_tmp(ioff,1,koff) -= tmp*f;
724  mat_tmp(ioff,0,koff) += tmp*f;
725  } else if (cell_id(ii,j-1,kk) < 0 || apy(ii,j-1,kk) == Real(0.0)) {
726  mat_tmp(ioff,0,koff) += tmp*(f+bflo);
727  } else {
728  mat_tmp(ioff, 0,koff) += tmp*(f+bf1);
729  mat_tmp(ioff,-1,koff) += tmp* bf2;
730  }
731  }
732  }
733 
734  cdir = Orientation(Direction::z, Orientation::low);
735  area = apz(i,j,k);
736 
737  if (area > Real(0.0))
738  {
739  int ioff, joff, ii, jj;
740  Real fracx, fracy;
741  Real f = fac[2];
742 
743  if (area != Real(1.0)) {
744  ioff = static_cast<int>(std::copysign(Real(1.0), fcz(i,j,k,0)));
745  joff = static_cast<int>(std::copysign(Real(1.0), fcz(i,j,k,1)));
746  ii = i+ioff;
747  jj = j+joff;
748  if (cell_id(ii,j,k-1) < 0 && cell_id(ii,j,k) < 0) {
749  fracx = Real(0.0);
750  } else {
751  fracx = std::abs(fcz(i,j,k,0));
752  }
753  if (cell_id(i,jj,k-1) < 0 && cell_id(i,jj,k) < 0) {
754  fracy = Real(0.0);
755  } else {
756  fracy = std::abs(fcz(i,j,k,1));
757  }
758  if (cell_id(ii,jj,k-1) < 0 && cell_id(ii,jj,k) < 0 && fracx*fracy > Real(0.0)) {
759  fracx = Real(0.0);
760  fracy = Real(0.0);
761  }
762  } else {
763  ioff = 0;
764  joff = 0;
765  ii = i;
766  jj = j;
767  fracx = Real(0.0);
768  fracy = Real(0.0);
769  }
770 
771  Real bf1, bf2, bflo;
772  detail::comp_bflo(bf1, bf2, bflo, sb, dx[2], bctype[cdir], bcl[cdir], bho);
773 
774  Real tmp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*area*b[2](i,j,k);
775  if (cell_id(i,j,k-1) >= 0) {
776  mat_tmp(0,0, 0) += tmp*f;
777  mat_tmp(0,0,-1) -= tmp*f;
778  } else if (cell_id(i,j,k+1) < 0 || apz(i,j,k+1) == Real(0.0)) {
779  mat_tmp(0,0,0) += tmp*(f+bflo);
780  } else {
781  mat_tmp(0,0,0) += tmp*(f+bf1);
782  mat_tmp(0,0,1) += tmp* bf2;
783  }
784 
785  if (fracx > Real(0.0)) {
786  tmp = fracx*(Real(1.0)-fracy)*area*b[2](ii,j,k);
787  if (cell_id(ii,j,k-1) >= 0 && cell_id(ii,j,k) >= 0) {
788  mat_tmp(ioff,0,-1) -= tmp*f;
789  mat_tmp(ioff,0, 0) += tmp*f;
790  } else if (cell_id(ii,j,k+1) < 0 || apz(ii,j,k+1) == Real(0.0)) {
791  mat_tmp(ioff,0,0) += tmp*(f+bflo);
792  } else {
793  mat_tmp(ioff,0,0) += tmp*(f+bf1);
794  mat_tmp(ioff,0,1) += tmp* bf2;
795  }
796  }
797 
798  if (fracy > Real(0.0)) {
799  tmp = fracy*(Real(1.0)-fracx)*area*b[2](i,jj,k);
800  if (cell_id(i,jj,k-1) >= 0 && cell_id(i,jj,k) >= 0) {
801  mat_tmp(0,joff,-1) -= tmp*f;
802  mat_tmp(0,joff, 0) += tmp*f;
803  } else if (cell_id(i,jj,k+1) < 0 || apz(i,jj,k+1) == Real(0.0)) {
804  mat_tmp(0,joff,0) += tmp*(f+bflo);
805  } else {
806  mat_tmp(0,joff,0) += tmp*(f+bf1);
807  mat_tmp(0,joff,1) += tmp* bf2;
808  }
809  }
810 
811  if (fracx > Real(0.0) && fracy > Real(0.0)) {
812  tmp = fracx*fracy*area*b[2](ii,jj,k);
813  if (cell_id(ii,jj,k-1) >= 0 && cell_id(ii,jj,k) >= 0) {
814  mat_tmp(ioff,joff,-1) -= tmp*f;
815  mat_tmp(ioff,joff, 0) += tmp*f;
816  } else if (cell_id(ii,jj,k+1) < 0 || apz(ii,jj,k+1) == Real(0.0)) {
817  mat_tmp(ioff,joff,0) += tmp*(f+bflo);
818  } else {
819  mat_tmp(ioff,joff,0) += tmp*(f+bf1);
820  mat_tmp(ioff,joff,1) += tmp* bf2;
821  }
822  }
823  }
824 
825  cdir = Orientation(Direction::z, Orientation::high);
826  area = apz(i,j,k+1);
827 
828  if (area > Real(0.0))
829  {
830  int ioff, joff, ii, jj;
831  Real fracx, fracy;
832  Real f = fac[2];
833 
834  if (area != Real(1.0)) {
835  ioff = static_cast<int>(std::copysign(Real(1.0), fcz(i,j,k+1,0)));
836  joff = static_cast<int>(std::copysign(Real(1.0), fcz(i,j,k+1,1)));
837  ii = i+ioff;
838  jj = j+joff;
839  if (cell_id(ii,j,k) < 0 && cell_id(ii,j,k+1) < 0) {
840  fracx = Real(0.0);
841  } else {
842  fracx = std::abs(fcz(i,j,k+1,0));
843  }
844  if (cell_id(i,jj,k) < 0 && cell_id(i,jj,k+1) < 0) {
845  fracy = Real(0.0);
846  } else {
847  fracy = std::abs(fcz(i,j,k+1,1));
848  }
849  if (cell_id(ii,jj,k) < 0 && cell_id(ii,jj,k+1) < 0 && fracx*fracy > Real(0.0)) {
850  fracx = Real(0.0);
851  fracy = Real(0.0);
852  }
853  } else {
854  ioff = 0;
855  joff = 0;
856  ii = i;
857  jj = j;
858  fracx = Real(0.0);
859  fracy = Real(0.0);
860  }
861 
862  Real bf1, bf2, bflo;
863  detail::comp_bflo(bf1, bf2, bflo, sb, dx[2], bctype[cdir], bcl[cdir], bho);
864 
865  Real tmp = (Real(1.0)-fracx)*(Real(1.0)-fracy)*area*b[2](i,j,k+1);
866  if (cell_id(i,j,k+1) >= 0) {
867  mat_tmp(0,0,0) += tmp*f;
868  mat_tmp(0,0,1) -= tmp*f;
869  } else if (cell_id(i,j,k-1) < 0 || apz(i,j,k-1) == Real(0.0)) {
870  mat_tmp(0,0,0) += tmp*(f+bflo);
871  } else {
872  mat_tmp(0,0, 0) += tmp*(f+bf1);
873  mat_tmp(0,0,-1) += tmp* bf2;
874  }
875 
876  if (fracx > Real(0.0)) {
877  tmp = fracx*(Real(1.0)-fracy)*area*b[2](ii,j,k+1);
878  if (cell_id(ii,j,k) >= 0 && cell_id(ii,j,k+1) >= 0) {
879  mat_tmp(ioff,0,0) += tmp*f;
880  mat_tmp(ioff,0,1) -= tmp*f;
881  } else if (cell_id(ii,j,k-1) < 0 || apz(ii,j,k-1) == Real(0.0)) {
882  mat_tmp(ioff,0,0) += tmp*(f+bflo);
883  } else {
884  mat_tmp(ioff,0, 0) += tmp*(f+bf1);
885  mat_tmp(ioff,0,-1) += tmp* bf2;
886  }
887  }
888 
889  if (fracy > Real(0.0)) {
890  tmp = fracy*(Real(1.0)-fracx)*area*b[2](i,jj,k+1);
891  if (cell_id(i,jj,k) >= 0 && cell_id(i,jj,k+1) >= 0) {
892  mat_tmp(0,joff,0) += tmp*f;
893  mat_tmp(0,joff,1) -= tmp*f;
894  } else if (cell_id(i,jj,k-1) < 0 || apz(i,jj,k-1) == Real(0.0)) {
895  mat_tmp(0,joff,0) += tmp*(f+bflo);
896  } else {
897  mat_tmp(0,joff, 0) += tmp*(f+bf1);
898  mat_tmp(0,joff,-1) += tmp* bf2;
899  }
900  }
901 
902  if (fracx > Real(0.0) && fracy > Real(0.0)) {
903  tmp = fracx*fracy*area*b[2](ii,jj,k+1);
904  if (cell_id(ii,jj,k+1) >= 0 && cell_id(ii,jj,k) >= 0) {
905  mat_tmp(ioff,joff,1) -= tmp*f;
906  mat_tmp(ioff,joff,0) += tmp*f;
907  } else if (cell_id(ii,jj,k-1) < 0 || apz(ii,jj,k-1) == Real(0.0)) {
908  mat_tmp(ioff,joff,0) += tmp*(f+bflo);
909  } else {
910  mat_tmp(ioff,joff, 0) += tmp*(f+bf1);
911  mat_tmp(ioff,joff,-1) += tmp* bf2;
912  }
913  }
914  }
915 
916  if (beb) {
917  Real anorm = std::sqrt((apx(i,j,k) - apx(i+1,j,k))*(apx(i,j,k) - apx(i+1,j,k))
918  + (apy(i,j,k) - apy(i,j+1,k))*(apy(i,j,k) - apy(i,j+1,k))
919  + (apz(i,j,k) - apz(i,j,k+1))*(apz(i,j,k) - apz(i,j,k+1)));
920 
921  Real anorminv = Real(1.0)/anorm;
922  Real anrmx = (apx(i,j,k) - apx(i+1,j,k))*anorminv;
923  Real anrmy = (apy(i,j,k) - apy(i,j+1,k))*anorminv;
924  Real anrmz = (apz(i,j,k) - apz(i,j,k+1))*anorminv;
925  Real sx = std::copysign(Real(1.0),anrmx);
926  Real sy = std::copysign(Real(1.0),anrmy);
927  Real sz = std::copysign(Real(1.0),anrmz);
928  Real bctx = bcen(i,j,k,0);
929  Real bcty = bcen(i,j,k,1);
930  Real bctz = bcen(i,j,k,2);
931  Real dg = get_dx_eb(vfrc(i,j,k)) / amrex::max(std::abs(anrmx), std::abs(anrmy), std::abs(anrmz));
932  Real gx = bctx - dg*anrmx;
933  Real gy = bcty - dg*anrmy;
934  Real gz = bctz - dg*anrmz;
935  int ioff = -static_cast<int>(sx);
936  int joff = -static_cast<int>(sy);
937  int koff = -static_cast<int>(sz);
938  gx *= sx;
939  gy *= sy;
940  gz *= sz;
941  Real gxy = gx*gy;
942  Real gxz = gx*gz;
943  Real gyz = gy*gz;
944  Real gxyz = gx*gy*gz;
945  Array1D<Real,0,7> phig1{Real(1.0) + gx + gy + gz + gxy + gxz + gyz + gxyz,
946  - gx - gxy - gxz - gxyz,
947  - gy - gxy - gyz - gxyz,
948  - gz - gxz - gyz - gxyz,
949  + gxy + gxyz,
950  + gxz + gxyz,
951  + gyz + gxyz,
952  - gxyz};
953  Array1D<Real,0,7> feb;
954  for (int ii=0; ii<8; ii++) {
955  feb(ii) = -phig1(ii) * (ba(i,j,k) * beb(i,j,k) / dg);
956  }
957  mat_tmp(0 , 0 , 0 ) -= feb(0)*fac[0];
958  mat_tmp(ioff, 0 , 0 ) -= feb(1)*fac[0];
959  mat_tmp(0 ,joff, 0 ) -= feb(2)*fac[0];
960  mat_tmp(0 , 0 ,koff) -= feb(3)*fac[0];
961  mat_tmp(ioff,joff, 0 ) -= feb(4)*fac[0];
962  mat_tmp(ioff, 0 ,koff) -= feb(5)*fac[0];
963  mat_tmp(0 ,joff,koff) -= feb(6)*fac[0];
964  mat_tmp(ioff,joff,koff) -= feb(7)*fac[0];
965  }
966 
967  for (int kk=-1; kk<=1; kk++) {
968  for (int jj=-1; jj<=1; jj++) {
969  for (int ii=-1; ii<=1; ii++) {
970  mat_tmp(ii,jj,kk) *= Real(1.0)/vfrc(i,j,k);
971  }}}
972 
973  if (sa != Real(0.0)) {
974  mat_tmp(0,0,0) += sa*a(i,j,k);
975  }
976  }
977  else
978  {
979  for (int kk=-1; kk<=1; kk++) {
980  for (int jj=-1; jj<=1; jj++) {
981  for (int ii=-1; ii<=1; ii++) {
982  mat_tmp(ii,jj,kk) = Real(0.0);
983  }}}
984  mat_tmp(0,0,0) = Real(1.0);
985  }
986 
987  diaginv(i,j,k) = Real(1.0)/mat_tmp(0,0,0);
988  for (int kk=-1; kk<=1; kk++) {
989  for (int jj=-1; jj<=1; jj++) {
990  for (int ii=-1; ii<=1; ii++) {
991  mat_tmp(ii,jj,kk) *= diaginv(i,j,k);
992  }}}
993  mat_tmp(0,0,0) = Real(1.0);
994 
995  ncols(i,j,k) = 0;
996  for (int m = 0; m < 27; ++m) {
997  ncols(i,j,k) += (sten[m] != Real(0.0));
998  }
999 }
1000 
1001 #endif
1002 
1003 }
1004 #endif
AMReX_EBFabFactory.H
AMReX_EBMultiFabUtil.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
AMReX_MLLinOp_K.H
AMReX_MultiCutFab.H
amrex::Orientation::low
@ low
Definition: AMReX_Orientation.H:34
amrex::Orientation::high
@ high
Definition: AMReX_Orientation.H:34
amrex::SundialsUserFun::f
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
amrex::detail::comp_bf
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void comp_bf(Real &bf1, Real &bf2, Real sb, Real h, int bct, Real bcl, int bho)
Definition: AMReX_Habec_K.H:11
amrex::detail::comp_bflo
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void comp_bflo(Real &bf1, Real &bf2, Real &bflo, Real sb, Real h, int bct, Real bcl, int bho)
Definition: AMReX_Habec_K.H:34
amrex
Definition: AMReX_Amr.cpp:49
amrex::habec_ijmat
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void habec_ijmat(GpuArray< Real, 2 *AMREX_SPACEDIM+1 > &sten, Array4< Int > const &ncols, Array4< Real > const &diaginv, int i, int j, int k, Array4< Int const > const &cell_id, Real sa, Array4< Real const > const &a, Real sb, GpuArray< Real, AMREX_SPACEDIM > const &dx, GpuArray< Array4< Real const >, AMREX_SPACEDIM > const &b, GpuArray< int, AMREX_SPACEDIM *2 > const &bctype, GpuArray< Real, AMREX_SPACEDIM *2 > const &bcl, int bho, Array4< int const > const &osm)
Definition: AMReX_Habec_2D_K.H:91
amrex::max
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
amrex::abs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
amrex::habec_cols
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void habec_cols(GpuArray< Int, 2 *AMREX_SPACEDIM+1 > &sten, int i, int j, int, Array4< Int const > const &cell_id)
Definition: AMReX_Habec_2D_K.H:171
amrex::habec_mat
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void habec_mat(GpuArray< Real, 2 *AMREX_SPACEDIM+1 > &sten, int i, int j, int k, Dim3 const &boxlo, Dim3 const &boxhi, Real sa, Array4< Real const > const &a, Real sb, GpuArray< Real, AMREX_SPACEDIM > const &dx, GpuArray< Array4< Real const >, AMREX_SPACEDIM > const &b, GpuArray< int, AMREX_SPACEDIM *2 > const &bctype, GpuArray< Real, AMREX_SPACEDIM *2 > const &bcl, int bho, GpuArray< Array4< int const >, AMREX_SPACEDIM *2 > const &msk, Array4< Real > const &diaginv)
Definition: AMReX_Habec_2D_K.H:15
amrex::sqrt
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373