Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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
10namespace amrex {
11
12template <typename T>
14void mllinop_apply_bc_x (int side, Box const& box, int blen,
15 Array4<T> const& phi,
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) {
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 }
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
74template <typename T>
76void mllinop_apply_bc_x (int side, int i, int j, int k, int blen,
77 Array4<T> const& phi,
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) {
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 }
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
117template <typename T>
119void 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 }
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
179template <typename T>
181void 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 }
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
222template <typename T>
224void 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 }
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
284template <typename T>
286void 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 }
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
327template <typename T>
329void 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 }
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
375template <typename T>
377void 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 }
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
408template <typename T>
410void 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 }
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
456template <typename T>
458void 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 }
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
489template <typename T>
491void 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 }
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
537template <typename T>
539void 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 }
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
573void 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 }
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
640void 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 }
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
692void 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 }
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
759void 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 }
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
811void 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 }
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
878void 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 }
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
930template <typename T>
932void 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
946template <typename T>
948void 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
962template <typename T>
964void 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
978template <typename T>
980void 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
993template <typename T>
995void 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
1009template <typename T>
1011void 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
1024template <typename T>
1026void 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
1040template <typename T>
1042void 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
1059Real 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
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 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 AMREX_FORCE_INLINE constexpr 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 constexpr 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_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