Block-Structured AMR Software Framework
AMReX_MLLinOp_K.H
Go to the documentation of this file.
1 #ifndef AMREX_MLLINOP_K_H_
2 #define AMREX_MLLINOP_K_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_FArrayBox.H>
6 #include <AMReX_BoundCond.H>
7 #include <AMReX_LO_BCTYPES.H>
8 #include <AMReX_LOUtil_K.H>
9 
10 namespace amrex {
11 
12 template <typename T>
14 void mllinop_apply_bc_x (int side, Box const& box, int blen,
15  Array4<T> const& phi,
16  Array4<int const> const& mask,
17  BoundCond bct, T bcl,
18  Array4<T const> const& bcval,
19  int maxorder, T dxinv, int inhomog, int icomp) noexcept
20 {
21  const auto lo = amrex::lbound(box);
22  const auto hi = amrex::ubound(box);
23  const int i = lo.x; // boundary cell
24  const int s = 1-2*side; // +1 for lo and -1 for hi
25  switch (bct) {
26  case AMREX_LO_NEUMANN:
27  {
28  for (int k = lo.z; k <= hi.z; ++k) {
29  for (int j = lo.y; j <= hi.y; ++j) {
30  if (mask(i,j,k) > 0) {
31  phi(i,j,k,icomp) = phi(i+s,j,k,icomp);
32  }
33  }
34  }
35  break;
36  }
38  {
39  for (int k = lo.z; k <= hi.z; ++k) {
40  for (int j = lo.y; j <= hi.y; ++j) {
41  if (mask(i,j,k) > 0) {
42  phi(i,j,k,icomp) = -phi(i+s,j,k,icomp);
43  }
44  }
45  }
46  break;
47  }
48  case AMREX_LO_DIRICHLET:
49  {
50  const int NX = amrex::min(blen+1, maxorder);
51  GpuArray<T,4> x{{-bcl * dxinv, T(0.5), T(1.5), T(2.5)}};
52  GpuArray<T,4> coef{};
53  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
54  for (int k = lo.z; k <= hi.z; ++k) {
55  for (int j = lo.y; j <= hi.y; ++j) {
56  if (mask(i,j,k) > 0) {
57  T tmp = T(0.0);
58  for (int m = 1; m < NX; ++m) {
59  tmp += phi(i+m*s,j,k,icomp) * coef[m];
60  }
61  phi(i,j,k,icomp) = tmp;
62  if (inhomog) {
63  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef[0];
64  }
65  }
66  }
67  }
68  break;
69  }
70  default: {}
71  }
72 }
73 
74 template <typename T>
76 void mllinop_apply_bc_x (int side, int i, int j, int k, int blen,
77  Array4<T> const& phi,
78  Array4<int const> const& mask,
79  BoundCond bct, T bcl,
80  Array4<T const> const& bcval,
81  int maxorder, T dxinv, int inhomog, int icomp) noexcept
82 {
83  if (mask(i,j,k) > 0) {
84  const int s = 1-2*side; // +1 for lo and -1 for hi
85  switch (bct) {
86  case AMREX_LO_NEUMANN:
87  {
88  phi(i,j,k,icomp) = phi(i+s,j,k,icomp);
89  break;
90  }
92  {
93  phi(i,j,k,icomp) = -phi(i+s,j,k,icomp);
94  break;
95  }
96  case AMREX_LO_DIRICHLET:
97  {
98  const int NX = amrex::min(blen+1, maxorder);
99  GpuArray<T,4> x{{-bcl * dxinv, T(0.5), T(1.5), T(2.5)}};
100  GpuArray<T,4> coef{};
101  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
102  T tmp = T(0.0);
103  for (int m = 1; m < NX; ++m) {
104  tmp += phi(i+m*s,j,k,icomp) * coef[m];
105  }
106  phi(i,j,k,icomp) = tmp;
107  if (inhomog) {
108  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef[0];
109  }
110  break;
111  }
112  default: {}
113  }
114  }
115 }
116 
117 template <typename T>
119 void mllinop_apply_bc_y (int side, Box const& box, int blen,
120  Array4<T> const& phi,
121  Array4<int const> const& mask,
122  BoundCond bct, T bcl,
123  Array4<T const> const& bcval,
124  int maxorder, T dyinv, int inhomog, int icomp) noexcept
125 {
126  const auto lo = amrex::lbound(box);
127  const auto hi = amrex::ubound(box);
128  const int j = lo.y; // boundary cell
129  const int s = 1-2*side; // +1 for lo and -1 for hi
130  switch (bct) {
131  case AMREX_LO_NEUMANN:
132  {
133  for (int k = lo.z; k <= hi.z; ++k) {
134  for (int i = lo.x; i <= hi.x; ++i) {
135  if (mask(i,j,k) > 0) {
136  phi(i,j,k,icomp) = phi(i,j+s,k,icomp);
137  }
138  }
139  }
140  break;
141  }
143  {
144  for (int k = lo.z; k <= hi.z; ++k) {
145  for (int i = lo.x; i <= hi.x; ++i) {
146  if (mask(i,j,k) > 0) {
147  phi(i,j,k,icomp) = -phi(i,j+s,k,icomp);
148  }
149  }
150  }
151  break;
152  }
153  case AMREX_LO_DIRICHLET:
154  {
155  const int NX = amrex::min(blen+1, maxorder);
156  GpuArray<T,4> x{{-bcl * dyinv, T(0.5), T(1.5), T(2.5)}};
157  GpuArray<T,4> coef{};
158  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
159  for (int k = lo.z; k <= hi.z; ++k) {
160  for (int i = lo.x; i <= hi.x; ++i) {
161  if (mask(i,j,k) > 0) {
162  T tmp = T(0.0);
163  for (int m = 1; m < NX; ++m) {
164  tmp += phi(i,j+m*s,k,icomp) * coef[m];
165  }
166  phi(i,j,k,icomp) = tmp;
167  if (inhomog) {
168  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef[0];
169  }
170  }
171  }
172  }
173  break;
174  }
175  default: {}
176  }
177 }
178 
179 template <typename T>
181 void mllinop_apply_bc_y (int side, int i, int j, int k, int blen,
182  Array4<T> const& phi,
183  Array4<int const> const& mask,
184  BoundCond bct, T bcl,
185  Array4<T const> const& bcval,
186  int maxorder, T dyinv, int inhomog, int icomp) noexcept
187 {
188  if (mask(i,j,k) > 0) {
189  const int s = 1-2*side; // +1 for lo and -1 for hi
190  switch (bct) {
191  case AMREX_LO_NEUMANN:
192  {
193  phi(i,j,k,icomp) = phi(i,j+s,k,icomp);
194  break;
195  }
197  {
198  phi(i,j,k,icomp) = -phi(i,j+s,k,icomp);
199  break;
200  }
201  case AMREX_LO_DIRICHLET:
202  {
203  const int NX = amrex::min(blen+1, maxorder);
204  GpuArray<T,4> x{{-bcl * dyinv, T(0.5), T(1.5), T(2.5)}};
205  GpuArray<T,4> coef{};
206  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
207  T tmp = T(0.0);
208  for (int m = 1; m < NX; ++m) {
209  tmp += phi(i,j+m*s,k,icomp) * coef[m];
210  }
211  phi(i,j,k,icomp) = tmp;
212  if (inhomog) {
213  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef[0];
214  }
215  break;
216  }
217  default: {}
218  }
219  }
220 }
221 
222 template <typename T>
224 void mllinop_apply_bc_z (int side, Box const& box, int blen,
225  Array4<T> const& phi,
226  Array4<int const> const& mask,
227  BoundCond bct, T bcl,
228  Array4<T const> const& bcval,
229  int maxorder, T dzinv, int inhomog, int icomp) noexcept
230 {
231  const auto lo = amrex::lbound(box);
232  const auto hi = amrex::ubound(box);
233  const int k = lo.z; // boundary cell
234  const int s = 1-2*side; // +1 for lo and -1 for hi
235  switch (bct) {
236  case AMREX_LO_NEUMANN:
237  {
238  for (int j = lo.y; j <= hi.y; ++j) {
239  for (int i = lo.x; i <= hi.x; ++i) {
240  if (mask(i,j,k) > 0) {
241  phi(i,j,k,icomp) = phi(i,j,k+s,icomp);
242  }
243  }
244  }
245  break;
246  }
248  {
249  for (int j = lo.y; j <= hi.y; ++j) {
250  for (int i = lo.x; i <= hi.x; ++i) {
251  if (mask(i,j,k) > 0) {
252  phi(i,j,k,icomp) = -phi(i,j,k+s,icomp);
253  }
254  }
255  }
256  break;
257  }
258  case AMREX_LO_DIRICHLET:
259  {
260  const int NX = amrex::min(blen+1, maxorder);
261  GpuArray<T,4> x{{-bcl * dzinv, T(0.5), T(1.5), T(2.5)}};
262  GpuArray<T,4> coef{};
263  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
264  for (int j = lo.y; j <= hi.y; ++j) {
265  for (int i = lo.x; i <= hi.x; ++i) {
266  if (mask(i,j,k) > 0) {
267  T tmp = T(0.0);
268  for (int m = 1; m < NX; ++m) {
269  tmp += phi(i,j,k+m*s,icomp) * coef[m];
270  }
271  phi(i,j,k,icomp) = tmp;
272  if (inhomog) {
273  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef[0];
274  }
275  }
276  }
277  }
278  break;
279  }
280  default: {}
281  }
282 }
283 
284 template <typename T>
286 void mllinop_apply_bc_z (int side, int i, int j, int k, int blen,
287  Array4<T> const& phi,
288  Array4<int const> const& mask,
289  BoundCond bct, T bcl,
290  Array4<T const> const& bcval,
291  int maxorder, T dzinv, int inhomog, int icomp) noexcept
292 {
293  if (mask(i,j,k) > 0) {
294  const int s = 1-2*side; // +1 for lo and -1 for hi
295  switch (bct) {
296  case AMREX_LO_NEUMANN:
297  {
298  phi(i,j,k,icomp) = phi(i,j,k+s,icomp);
299  break;
300  }
302  {
303  phi(i,j,k,icomp) = -phi(i,j,k+s,icomp);
304  break;
305  }
306  case AMREX_LO_DIRICHLET:
307  {
308  const int NX = amrex::min(blen+1, maxorder);
309  GpuArray<T,4> x{{-bcl * dzinv, T(0.5), T(1.5), T(2.5)}};
310  GpuArray<T,4> coef{};
311  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
312  T tmp = T(0.0);
313  for (int m = 1; m < NX; ++m) {
314  tmp += phi(i,j,k+m*s,icomp) * coef[m];
315  }
316  phi(i,j,k,icomp) = tmp;
317  if (inhomog) {
318  phi(i,j,k,icomp) += bcval(i,j,k,icomp)*coef[0];
319  }
320  break;
321  }
322  default: {}
323  }
324  }
325 }
326 
327 template <typename T>
329 void mllinop_comp_interp_coef0_x (int side, Box const& box, int blen,
330  Array4<T> const& f,
331  Array4<int const> const& mask,
332  BoundCond bct, T bcl,
333  int maxorder, T dxinv, int icomp) noexcept
334 {
335  const auto lo = amrex::lbound(box);
336  const auto hi = amrex::ubound(box);
337  const int ib = lo.x; // boundary cell
338  const int ii = lo.x + (1-2*side); // interior cell
339  switch (bct) {
340  case AMREX_LO_NEUMANN:
341  {
342  for (int k = lo.z; k <= hi.z; ++k) {
343  for (int j = lo.y; j <= hi.y; ++j) {
344  f(ii,j,k,icomp) = T(1.0);
345  }
346  }
347  break;
348  }
350  {
351  for (int k = lo.z; k <= hi.z; ++k) {
352  for (int j = lo.y; j <= hi.y; ++j) {
353  f(ii,j,k,icomp) = (mask(ib,j,k) > 0) ? T(1.0) : T(0.0);
354  }
355  }
356  break;
357  }
358  case AMREX_LO_DIRICHLET:
359  {
360  const int NX = amrex::min(blen+1, maxorder);
361  GpuArray<T,4> x{{-bcl * dxinv, T(0.5), T(1.5), T(2.5)}};
362  GpuArray<T,4> coef{};
363  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
364  for (int k = lo.z; k <= hi.z; ++k) {
365  for (int j = lo.y; j <= hi.y; ++j) {
366  f(ii,j,k,icomp) = (mask(ib,j,k) > 0) ? coef[1] : T(0.0);
367  }
368  }
369  break;
370  }
371  default: {}
372  }
373 }
374 
375 template <typename T>
377 void mllinop_comp_interp_coef0_x (int side, int i, int j, int k, int blen,
378  Array4<T> const& f,
379  Array4<int const> const& mask,
380  BoundCond bct, T bcl,
381  int maxorder, T dxinv, int icomp) noexcept
382 {
383  const int ii = i + (1-2*side); // interior cell
384  switch (bct) {
385  case AMREX_LO_NEUMANN:
386  {
387  f(ii,j,k,icomp) = T(1.0);
388  break;
389  }
391  {
392  f(ii,j,k,icomp) = (mask(i,j,k) > 0) ? T(1.0) : T(0.0);
393  break;
394  }
395  case AMREX_LO_DIRICHLET:
396  {
397  const int NX = amrex::min(blen+1, maxorder);
398  GpuArray<T,4> x{{-bcl * dxinv, T(0.5), T(1.5), T(2.5)}};
399  GpuArray<T,4> coef{};
400  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
401  f(ii,j,k,icomp) = (mask(i,j,k) > 0) ? coef[1] : T(0.0);
402  break;
403  }
404  default: {}
405  }
406 }
407 
408 template <typename T>
410 void mllinop_comp_interp_coef0_y (int side, Box const& box, int blen,
411  Array4<T> const& f,
412  Array4<int const> const& mask,
413  BoundCond bct, T bcl,
414  int maxorder, T dyinv, int icomp) noexcept
415 {
416  const auto lo = amrex::lbound(box);
417  const auto hi = amrex::ubound(box);
418  const int jb = lo.y; // boundary cell
419  const int ji = lo.y + (1-2*side); // interior cell
420  switch (bct) {
421  case AMREX_LO_NEUMANN:
422  {
423  for (int k = lo.z; k <= hi.z; ++k) {
424  for (int i = lo.x; i <= hi.x; ++i) {
425  f(i,ji,k,icomp) = T(1.0);
426  }
427  }
428  break;
429  }
431  {
432  for (int k = lo.z; k <= hi.z; ++k) {
433  for (int i = lo.x; i <= hi.x; ++i) {
434  f(i,ji,k,icomp) = (mask(i,jb,k) > 0) ? T(1.0) : T(0.0);
435  }
436  }
437  break;
438  }
439  case AMREX_LO_DIRICHLET:
440  {
441  const int NX = amrex::min(blen+1, maxorder);
442  GpuArray<T,4> x{{-bcl * dyinv, T(0.5), T(1.5), T(2.5)}};
443  GpuArray<T,4> coef{};
444  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
445  for (int k = lo.z; k <= hi.z; ++k) {
446  for (int i = lo.x; i <= hi.x; ++i) {
447  f(i,ji,k,icomp) = (mask(i,jb,k) > 0) ? coef[1] : T(0.0);
448  }
449  }
450  break;
451  }
452  default: {}
453  }
454 }
455 
456 template <typename T>
458 void mllinop_comp_interp_coef0_y (int side, int i, int j, int k, int blen,
459  Array4<T> const& f,
460  Array4<int const> const& mask,
461  BoundCond bct, T bcl,
462  int maxorder, T dyinv, int icomp) noexcept
463 {
464  const int ji = j + (1-2*side); // interior cell
465  switch (bct) {
466  case AMREX_LO_NEUMANN:
467  {
468  f(i,ji,k,icomp) = T(1.0);
469  break;
470  }
472  {
473  f(i,ji,k,icomp) = (mask(i,j,k) > 0) ? T(1.0) : T(0.0);
474  break;
475  }
476  case AMREX_LO_DIRICHLET:
477  {
478  const int NX = amrex::min(blen+1, maxorder);
479  GpuArray<T,4> x{{-bcl * dyinv, T(0.5), T(1.5), T(2.5)}};
480  GpuArray<T,4> coef{};
481  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
482  f(i,ji,k,icomp) = (mask(i,j,k) > 0) ? coef[1] : T(0.0);
483  break;
484  }
485  default: {}
486  }
487 }
488 
489 template <typename T>
491 void mllinop_comp_interp_coef0_z (int side, Box const& box, int blen,
492  Array4<T> const& f,
493  Array4<int const> const& mask,
494  BoundCond bct, T bcl,
495  int maxorder, T dzinv, int icomp) noexcept
496 {
497  const auto lo = amrex::lbound(box);
498  const auto hi = amrex::ubound(box);
499  const int kb = lo.z; // bound cell
500  const int ki = lo.z + (1-2*side); // interior cell
501  switch (bct) {
502  case AMREX_LO_NEUMANN:
503  {
504  for (int j = lo.y; j <= hi.y; ++j) {
505  for (int i = lo.x; i <= hi.x; ++i) {
506  f(i,j,ki,icomp) = T(1.0);
507  }
508  }
509  break;
510  }
512  {
513  for (int j = lo.y; j <= hi.y; ++j) {
514  for (int i = lo.x; i <= hi.x; ++i) {
515  f(i,j,ki,icomp) = (mask(i,j,kb) > 0) ? T(1.0) : T(0.0);
516  }
517  }
518  break;
519  }
520  case AMREX_LO_DIRICHLET:
521  {
522  const int NX = amrex::min(blen+1, maxorder);
523  GpuArray<T,4> x{{-bcl * dzinv, T(0.5), T(1.5), T(2.5)}};
524  GpuArray<T,4> coef{};
525  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
526  for (int j = lo.y; j <= hi.y; ++j) {
527  for (int i = lo.x; i <= hi.x; ++i) {
528  f(i,j,ki,icomp) = (mask(i,j,kb) > 0) ? coef[1] : T(0.0);
529  }
530  }
531  break;
532  }
533  default: {}
534  }
535 }
536 
537 template <typename T>
539 void mllinop_comp_interp_coef0_z (int side, int i, int j, int k, int blen,
540  Array4<T> const& f,
541  Array4<int const> const& mask,
542  BoundCond bct, T bcl,
543  int maxorder, T dzinv, int icomp) noexcept
544 {
545  const int ki = k + (1-2*side); // interior cell
546  switch (bct) {
547  case AMREX_LO_NEUMANN:
548  {
549  f(i,j,ki,icomp) = T(1.0);
550  break;
551  }
553  {
554  f(i,j,ki,icomp) = (mask(i,j,k) > 0) ? T(1.0) : T(0.0);
555  break;
556  }
557  case AMREX_LO_DIRICHLET:
558  {
559  const int NX = amrex::min(blen+1, maxorder);
560  GpuArray<T,4> x{{-bcl * dzinv, T(0.5), T(1.5), T(2.5)}};
561  GpuArray<T,4> coef{};
562  poly_interp_coeff(-T(0.5), x.data(), NX, coef.data());
563  f(i,j,ki,icomp) = (mask(i,j,k) > 0) ? coef[1] : T(0.0);
564  break;
565  }
566  default: {}
567  }
568 }
569 
570 #ifdef AMREX_USE_EB
571 
573 void mllinop_comp_interp_coef0_x_eb (int side, Box const& box, int blen,
574  Array4<Real> const& f,
575  Array4<int const> const& mask,
576  Array4<Real const> const& area,
577  BoundCond bct, Real bcl,
578  int maxorder, Real dxinv, int icomp) noexcept
579 {
580  const auto lo = amrex::lbound(box);
581  const auto hi = amrex::ubound(box);
582  const int ib = lo.x; // boundary cell
583  const int s = 1-2*side; // +1 for lo and -1 for hi
584  const int ii = lo.x + s; // interior cell
585  switch (bct) {
586  case AMREX_LO_NEUMANN:
587  {
588  for (int k = lo.z; k <= hi.z; ++k) {
589  for (int j = lo.y; j <= hi.y; ++j) {
590  f(ii,j,k,icomp) = Real(1.0);
591  }
592  }
593  break;
594  }
596  {
597  for (int k = lo.z; k <= hi.z; ++k) {
598  for (int j = lo.y; j <= hi.y; ++j) {
599  f(ii,j,k,icomp) = (mask(ib,j,k) > 0) ? Real(1.0) : Real(0.0);
600  }
601  }
602  break;
603  }
604  case AMREX_LO_DIRICHLET:
605  {
606  const int NX = amrex::min(blen+1, maxorder);
607  GpuArray<Real,4> x{{-bcl * dxinv, Real(0.5), Real(1.5), Real(2.5)}};
608  Array2D<Real, 0, 3, 0, 2> coef{};
609  for (int r = 0; r <= maxorder-2; ++r) {
610  poly_interp_coeff(-Real(0.5), x.data(), r+2, &coef(0,r));
611  }
612  for (int k = lo.z; k <= hi.z; ++k) {
613  for (int j = lo.y; j <= hi.y; ++j) {
614  int order = 1;
615  if (mask(ib,j,k) > 0) {
616  bool has_cutfaces = false;
617  for (int r = 0; r <= NX-2; ++r) {
618  Real a = area(ii+side+s*r,j,k);
619  if (a > Real(0.0)) {
620  ++order;
621  if (a < Real(1.0)) {
622  has_cutfaces = true;
623  }
624  } else {
625  break;
626  }
627  }
628  if (has_cutfaces) { order = amrex::min(2,order); }
629  }
630  f(ii,j,k,icomp) = (order==1) ? Real(0.0) : coef(1,order-2);
631  }
632  }
633  break;
634  }
635  default: {}
636  }
637 }
638 
640 void mllinop_comp_interp_coef0_x_eb (int side, int i, int j, int k, int blen,
641  Array4<Real> const& f,
642  Array4<int const> const& mask,
643  Array4<Real const> const& area,
644  BoundCond bct, Real bcl,
645  int maxorder, Real dxinv, int icomp) noexcept
646 {
647  const int s = 1-2*side; // +1 for lo and -1 for hi
648  const int ii = i + s; // interior cell
649  switch (bct) {
650  case AMREX_LO_NEUMANN:
651  {
652  f(ii,j,k,icomp) = Real(1.0);
653  break;
654  }
656  {
657  f(ii,j,k,icomp) = (mask(i,j,k) > 0) ? Real(1.0) : Real(0.0);
658  break;
659  }
660  case AMREX_LO_DIRICHLET:
661  {
662  const int NX = amrex::min(blen+1, maxorder);
663  GpuArray<Real,4> x{{-bcl * dxinv, Real(0.5), Real(1.5), Real(2.5)}};
664  Array2D<Real, 0, 3, 0, 2> coef{};
665  for (int r = 0; r <= maxorder-2; ++r) {
666  poly_interp_coeff(-Real(0.5), x.data(), r+2, &coef(0,r));
667  }
668  int order = 1;
669  if (mask(i,j,k) > 0) {
670  bool has_cutfaces = false;
671  for (int r = 0; r <= NX-2; ++r) {
672  Real a = area(ii+side+s*r,j,k);
673  if (a > Real(0.0)) {
674  ++order;
675  if (a < Real(1.0)) {
676  has_cutfaces = true;
677  }
678  } else {
679  break;
680  }
681  }
682  if (has_cutfaces) { order = amrex::min(2,order); }
683  }
684  f(ii,j,k,icomp) = (order==1) ? Real(0.0) : coef(1,order-2);
685  break;
686  }
687  default: {}
688  }
689 }
690 
692 void mllinop_comp_interp_coef0_y_eb (int side, Box const& box, int blen,
693  Array4<Real> const& f,
694  Array4<int const> const& mask,
695  Array4<Real const> const& area,
696  BoundCond bct, Real bcl,
697  int maxorder, Real dyinv, int icomp) noexcept
698 {
699  const auto lo = amrex::lbound(box);
700  const auto hi = amrex::ubound(box);
701  const int jb = lo.y; // boundary cell
702  const int s = 1-2*side; // +1 for lo and -1 for hi
703  const int ji = lo.y + s; // interior cell
704  switch (bct) {
705  case AMREX_LO_NEUMANN:
706  {
707  for (int k = lo.z; k <= hi.z; ++k) {
708  for (int i = lo.x; i <= hi.x; ++i) {
709  f(i,ji,k,icomp) = Real(1.0);
710  }
711  }
712  break;
713  }
715  {
716  for (int k = lo.z; k <= hi.z; ++k) {
717  for (int i = lo.x; i <= hi.x; ++i) {
718  f(i,ji,k,icomp) = (mask(i,jb,k) > 0) ? Real(1.0) : Real(0.0);
719  }
720  }
721  break;
722  }
723  case AMREX_LO_DIRICHLET:
724  {
725  const int NX = amrex::min(blen+1, maxorder);
726  GpuArray<Real,4> x{{-bcl * dyinv, Real(0.5), Real(1.5), Real(2.5)}};
727  Array2D<Real, 0, 3, 0, 2> coef{};
728  for (int r = 0; r <= maxorder-2; ++r) {
729  poly_interp_coeff(-Real(0.5), x.data(), r+2, &coef(0,r));
730  }
731  for (int k = lo.z; k <= hi.z; ++k) {
732  for (int i = lo.x; i <= hi.x; ++i) {
733  int order = 1;
734  if (mask(i,jb,k) > 0) {
735  bool has_cutfaces = false;
736  for (int r = 0; r <= NX-2; ++r) {
737  Real a = area(i,ji+side+s*r,k);
738  if (a > Real(0.0)) {
739  ++order;
740  if (a < Real(1.0)) {
741  has_cutfaces = true;
742  }
743  } else {
744  break;
745  }
746  }
747  if (has_cutfaces) { order = amrex::min(2,order); }
748  }
749  f(i,ji,k,icomp) = (order==1) ? Real(0.0) : coef(1,order-2);
750  }
751  }
752  break;
753  }
754  default: {}
755  }
756 }
757 
759 void mllinop_comp_interp_coef0_y_eb (int side, int i, int j, int k, int blen,
760  Array4<Real> const& f,
761  Array4<int const> const& mask,
762  Array4<Real const> const& area,
763  BoundCond bct, Real bcl,
764  int maxorder, Real dyinv, int icomp) noexcept
765 {
766  const int s = 1-2*side; // +1 for lo and -1 for hi
767  const int ji = j + s; // interior cell
768  switch (bct) {
769  case AMREX_LO_NEUMANN:
770  {
771  f(i,ji,k,icomp) = Real(1.0);
772  break;
773  }
775  {
776  f(i,ji,k,icomp) = (mask(i,j,k) > 0) ? Real(1.0) : Real(0.0);
777  break;
778  }
779  case AMREX_LO_DIRICHLET:
780  {
781  const int NX = amrex::min(blen+1, maxorder);
782  GpuArray<Real,4> x{{-bcl * dyinv, Real(0.5), Real(1.5), Real(2.5)}};
783  Array2D<Real, 0, 3, 0, 2> coef{};
784  for (int r = 0; r <= maxorder-2; ++r) {
785  poly_interp_coeff(-Real(0.5), x.data(), r+2, &coef(0,r));
786  }
787  int order = 1;
788  if (mask(i,j,k) > 0) {
789  bool has_cutfaces = false;
790  for (int r = 0; r <= NX-2; ++r) {
791  Real a = area(i,ji+side+s*r,k);
792  if (a > Real(0.0)) {
793  ++order;
794  if (a < Real(1.0)) {
795  has_cutfaces = true;
796  }
797  } else {
798  break;
799  }
800  }
801  if (has_cutfaces) { order = amrex::min(2,order); }
802  }
803  f(i,ji,k,icomp) = (order==1) ? Real(0.0) : coef(1,order-2);
804  break;
805  }
806  default: {}
807  }
808 }
809 
811 void mllinop_comp_interp_coef0_z_eb (int side, Box const& box, int blen,
812  Array4<Real> const& f,
813  Array4<int const> const& mask,
814  Array4<Real const> const& area,
815  BoundCond bct, Real bcl,
816  int maxorder, Real dzinv, int icomp) noexcept
817 {
818  const auto lo = amrex::lbound(box);
819  const auto hi = amrex::ubound(box);
820  const int kb = lo.z; // bound cell
821  const int s = 1-2*side; // +1 for lo and -1 for hi
822  const int ki = lo.z + s; // interior cell
823  switch (bct) {
824  case AMREX_LO_NEUMANN:
825  {
826  for (int j = lo.y; j <= hi.y; ++j) {
827  for (int i = lo.x; i <= hi.x; ++i) {
828  f(i,j,ki,icomp) = Real(1.0);
829  }
830  }
831  break;
832  }
834  {
835  for (int j = lo.y; j <= hi.y; ++j) {
836  for (int i = lo.x; i <= hi.x; ++i) {
837  f(i,j,ki,icomp) = (mask(i,j,kb) > 0) ? Real(1.0) : Real(0.0);
838  }
839  }
840  break;
841  }
842  case AMREX_LO_DIRICHLET:
843  {
844  const int NX = amrex::min(blen+1, maxorder);
845  GpuArray<Real,4> x{{-bcl * dzinv, Real(0.5), Real(1.5), Real(2.5)}};
846  Array2D<Real, 0, 3, 0, 2> coef{};
847  for (int r = 0; r <= maxorder-2; ++r) {
848  poly_interp_coeff(-Real(0.5), x.data(), r+2, &coef(0,r));
849  }
850  for (int j = lo.y; j <= hi.y; ++j) {
851  for (int i = lo.x; i <= hi.x; ++i) {
852  int order = 1;
853  if (mask(i,j,kb) > 0) {
854  bool has_cutfaces = false;
855  for (int r = 0; r <= NX-2; ++r) {
856  Real a = area(i,j,ki+side+s*r);
857  if (a > Real(0.0)) {
858  ++order;
859  if (a < Real(1.0)) {
860  has_cutfaces = true;
861  }
862  } else {
863  break;
864  }
865  }
866  if (has_cutfaces) { order = amrex::min(2,order); }
867  }
868  f(i,j,ki,icomp) = (order==1) ? Real(0.0) : coef(1,order-2);
869  }
870  }
871  break;
872  }
873  default: {}
874  }
875 }
876 
878 void mllinop_comp_interp_coef0_z_eb (int side, int i, int j, int k, int blen,
879  Array4<Real> const& f,
880  Array4<int const> const& mask,
881  Array4<Real const> const& area,
882  BoundCond bct, Real bcl,
883  int maxorder, Real dzinv, int icomp) noexcept
884 {
885  const int s = 1-2*side; // +1 for lo and -1 for hi
886  const int ki = k + s; // interior cell
887  switch (bct) {
888  case AMREX_LO_NEUMANN:
889  {
890  f(i,j,ki,icomp) = Real(1.0);
891  break;
892  }
894  {
895  f(i,j,ki,icomp) = (mask(i,j,k) > 0) ? Real(1.0) : Real(0.0);
896  break;
897  }
898  case AMREX_LO_DIRICHLET:
899  {
900  const int NX = amrex::min(blen+1, maxorder);
901  GpuArray<Real,4> x{{-bcl * dzinv, Real(0.5), Real(1.5), Real(2.5)}};
902  Array2D<Real, 0, 3, 0, 2> coef{};
903  for (int r = 0; r <= maxorder-2; ++r) {
904  poly_interp_coeff(-Real(0.5), x.data(), r+2, &coef(0,r));
905  }
906  int order = 1;
907  if (mask(i,j,k) > 0) {
908  bool has_cutfaces = false;
909  for (int r = 0; r <= NX-2; ++r) {
910  Real a = area(i,j,ki+side+s*r);
911  if (a > Real(0.0)) {
912  ++order;
913  if (a < Real(1.0)) {
914  has_cutfaces = true;
915  }
916  } else {
917  break;
918  }
919  }
920  if (has_cutfaces) { order = amrex::min(2,order); }
921  }
922  f(i,j,ki,icomp) = (order==1) ? Real(0.0) : coef(1,order-2);
923  break;
924  }
925  default: {}
926  }
927 }
928 #endif
929 
930 template <typename T>
932 void mllinop_apply_innu_xlo (int i, int j, int k,
933  Array4<T> const& rhs,
934  Array4<int const> const& mask,
935  Array4<T const> const& bcoef,
936  BoundCond bct, T /*bcl*/,
937  Array4<T const> const& bcval,
938  T fac, bool has_bcoef, int icomp) noexcept
939 {
940  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
941  T b = (has_bcoef) ? bcoef(i+1,j,k,icomp) : T(1.0);
942  rhs(i+1,j,k,icomp) -= fac*b*bcval(i,j,k,icomp);
943  }
944 }
945 
946 template <typename T>
948 void mllinop_apply_innu_xhi (int i, int j, int k,
949  Array4<T> const& rhs,
950  Array4<int const> const& mask,
951  Array4<T const> const& bcoef,
952  BoundCond bct, T /*bcl*/,
953  Array4<T const> const& bcval,
954  T fac, bool has_bcoef, int icomp) noexcept
955 {
956  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
957  T b = (has_bcoef) ? bcoef(i,j,k,icomp) : T(1.0);
958  rhs(i-1,j,k,icomp) += fac*b*bcval(i,j,k,icomp);
959  }
960 }
961 
962 template <typename T>
964 void mllinop_apply_innu_ylo (int i, int j, int k,
965  Array4<T> const& rhs,
966  Array4<int const> const& mask,
967  Array4<T const> const& bcoef,
968  BoundCond bct, T /*bcl*/,
969  Array4<T const> const& bcval,
970  T fac, bool has_bcoef, int icomp) noexcept
971 {
972  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
973  T b = (has_bcoef) ? bcoef(i,j+1,k,icomp) : T(1.0);
974  rhs(i,j+1,k,icomp) -= fac*b*bcval(i,j,k,icomp);
975  }
976 }
977 
978 template <typename T>
980 void mllinop_apply_innu_ylo_m (int i, int j, int k,
981  Array4<T> const& rhs,
982  Array4<int const> const& mask,
983  BoundCond bct, T /*bcl*/,
984  Array4<T const> const& bcval,
985  T fac, T xlo, T dx, int icomp) noexcept
986 {
987  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
988  T b = xlo + (i+T(0.5))*dx;
989  rhs(i,j+1,k,icomp) -= fac*b*bcval(i,j,k,icomp);
990  }
991 }
992 
993 template <typename T>
995 void mllinop_apply_innu_yhi (int i, int j, int k,
996  Array4<T> const& rhs,
997  Array4<int const> const& mask,
998  Array4<T const> const& bcoef,
999  BoundCond bct, T /*bcl*/,
1000  Array4<T const> const& bcval,
1001  T fac, bool has_bcoef, int icomp) noexcept
1002 {
1003  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
1004  T b = (has_bcoef) ? bcoef(i,j,k,icomp) : T(1.0);
1005  rhs(i,j-1,k,icomp) += fac*b*bcval(i,j,k,icomp);
1006  }
1007 }
1008 
1009 template <typename T>
1011 void mllinop_apply_innu_yhi_m (int i, int j, int k,
1012  Array4<T> const& rhs,
1013  Array4<int const> const& mask,
1014  BoundCond bct, T /*bcl*/,
1015  Array4<T const> const& bcval,
1016  T fac, T xlo, T dx, int icomp) noexcept
1017 {
1018  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
1019  T b = xlo + (i+T(0.5))*dx;
1020  rhs(i,j-1,k,icomp) += fac*b*bcval(i,j,k,icomp);
1021  }
1022 }
1023 
1024 template <typename T>
1026 void mllinop_apply_innu_zlo (int i, int j, int k,
1027  Array4<T> const& rhs,
1028  Array4<int const> const& mask,
1029  Array4<T const> const& bcoef,
1030  BoundCond bct, T /*bcl*/,
1031  Array4<T const> const& bcval,
1032  T fac, bool has_bcoef, int icomp) noexcept
1033 {
1034  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
1035  T b = (has_bcoef) ? bcoef(i,j,k+1,icomp) : T(1.0);
1036  rhs(i,j,k+1,icomp) -= fac*b*bcval(i,j,k,icomp);
1037  }
1038 }
1039 
1040 template <typename T>
1042 void mllinop_apply_innu_zhi (int i, int j, int k,
1043  Array4<T> const& rhs,
1044  Array4<int const> const& mask,
1045  Array4<T const> const& bcoef,
1046  BoundCond bct, T /*bcl*/,
1047  Array4<T const> const& bcval,
1048  T fac, bool has_bcoef, int icomp) noexcept
1049 {
1050  if (bct == AMREX_LO_NEUMANN && mask(i,j,k) == 2) {
1051  T b = (has_bcoef) ? bcoef(i,j,k,icomp) : T(1.0);
1052  rhs(i,j,k-1,icomp) += fac*b*bcval(i,j,k,icomp);
1053  }
1054 }
1055 
1056 #ifdef AMREX_USE_EB
1057 
1059 Real get_dx_eb (Real kappa) noexcept {
1060  return amrex::max(Real(0.3),(kappa*kappa-Real(0.25))/(Real(2.0)*kappa));
1061 }
1062 
1063 #endif
1064 
1065 }
1066 
1067 #endif
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
Array4< int const > mask
Definition: AMReX_InterpFaceRegister.cpp:93
#define AMREX_LO_NEUMANN
Definition: AMReX_LO_BCTYPES.H:6
#define AMREX_LO_DIRICHLET
Definition: AMReX_LO_BCTYPES.H:5
#define AMREX_LO_REFLECT_ODD
Definition: AMReX_LO_BCTYPES.H:7
Maintain an identifier for boundary condition types.
Definition: AMReX_BoundCond.H:20
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
Definition: AMReX_Amr.cpp:49
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_zlo(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, Array4< T const > const &bcoef, BoundCond bct, T, Array4< T const > const &bcval, T fac, bool has_bcoef, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:1026
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_xhi(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, Array4< T const > const &bcoef, BoundCond bct, T, Array4< T const > const &bcval, T fac, bool has_bcoef, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:948
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_ylo(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, Array4< T const > const &bcoef, BoundCond bct, T, Array4< T const > const &bcval, T fac, bool has_bcoef, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:964
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_zhi(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, Array4< T const > const &bcoef, BoundCond bct, T, Array4< T const > const &bcval, T fac, bool has_bcoef, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:1042
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void poly_interp_coeff(T xInt, T const *AMREX_RESTRICT x, int N, T *AMREX_RESTRICT c) noexcept
Definition: AMReX_LOUtil_K.H:24
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_xlo(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, Array4< T const > const &bcoef, BoundCond bct, T, Array4< T const > const &bcval, T fac, bool has_bcoef, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:932
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
Definition: AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_ylo_m(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, BoundCond bct, T, Array4< T const > const &bcval, T fac, T xlo, T dx, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:980
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_yhi(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, Array4< T const > const &bcoef, BoundCond bct, T, Array4< T const > const &bcval, T fac, bool has_bcoef, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:995
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_innu_yhi_m(int i, int j, int k, Array4< T > const &rhs, Array4< int const > const &mask, BoundCond bct, T, Array4< T const > const &bcval, T fac, T xlo, T dx, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:1011
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_bc_x(int side, Box const &box, int blen, Array4< T > const &phi, Array4< int const > const &mask, BoundCond bct, T bcl, Array4< T const > const &bcval, int maxorder, T dxinv, int inhomog, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:14
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_comp_interp_coef0_x(int side, Box const &box, int blen, Array4< T > const &f, Array4< int const > const &mask, BoundCond bct, T bcl, int maxorder, T dxinv, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:329
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_bc_z(int side, Box const &box, int blen, Array4< T > const &phi, Array4< int const > const &mask, BoundCond bct, T bcl, Array4< T const > const &bcval, int maxorder, T dzinv, int inhomog, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:224
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_comp_interp_coef0_y(int side, Box const &box, int blen, Array4< T > const &f, Array4< int const > const &mask, BoundCond bct, T bcl, int maxorder, T dyinv, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:410
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_comp_interp_coef0_z(int side, Box const &box, int blen, Array4< T > const &f, Array4< int const > const &mask, BoundCond bct, T bcl, int maxorder, T dzinv, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:491
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mllinop_apply_bc_y(int side, Box const &box, int blen, Array4< T > const &phi, Array4< int const > const &mask, BoundCond bct, T bcl, Array4< T const > const &bcval, int maxorder, T dyinv, int inhomog, int icomp) noexcept
Definition: AMReX_MLLinOp_K.H:119
Definition: AMReX_Array4.H:61
Definition: AMReX_Array.H:34