1 #ifndef AMREX_MLNODETENSORLAP_3D_K_H_
2 #define AMREX_MLNODETENSORLAP_3D_K_H_
3 #include <AMReX_Config.H>
7 namespace mlndts_detail {
11 int ic,
int jc,
int kc) noexcept
13 return (
crse(ic,jc,kc)+
crse(ic+1,jc,kc))*Real(0.5);
18 int ic,
int jc,
int kc) noexcept
20 return (
crse(ic,jc,kc)+
crse(ic,jc+1,kc))*Real(0.5);
25 int ic,
int jc,
int kc) noexcept
27 return (
crse(ic,jc,kc)+
crse(ic,jc,kc+1))*Real(0.5);
32 int ic,
int jc,
int kc) noexcept
42 int ic,
int jc,
int kc) noexcept
52 int ic,
int jc,
int kc) noexcept
63 Array4<Real const>
const&
crse, Array4<int const>
const& msk) noexcept
65 using namespace mlndts_detail;
71 bool i_is_odd = (ic*2 != i);
72 bool j_is_odd = (jc*2 != j);
73 bool k_is_odd = (kc*2 != k);
74 if (i_is_odd && j_is_odd && k_is_odd) {
82 }
else if (j_is_odd && k_is_odd) {
85 }
else if (i_is_odd && k_is_odd) {
88 }
else if (i_is_odd && j_is_odd) {
91 }
else if (i_is_odd) {
94 }
else if (j_is_odd) {
97 }
else if (k_is_odd) {
109 Array4<Real const>
const&
crse, Array4<int const>
const& msk,
110 int semi_dir) noexcept
112 using namespace mlndts_detail;
118 bool j_is_odd = (jc*2 != j);
119 bool k_is_odd = (kc*2 != k);
120 if (j_is_odd && k_is_odd) {
123 }
else if (j_is_odd) {
126 }
else if (k_is_odd) {
133 }
else if (semi_dir == 1) {
136 bool i_is_odd = (ic*2 != i);
137 bool k_is_odd = (kc*2 != k);
138 if (i_is_odd && k_is_odd) {
141 }
else if (i_is_odd) {
144 }
else if (k_is_odd) {
154 bool i_is_odd = (ic*2 != i);
155 bool j_is_odd = (jc*2 != j);
156 if (i_is_odd && j_is_odd) {
159 }
else if (i_is_odd) {
162 }
else if (j_is_odd) {
178 y(i,j,k) = Real(0.0);
180 y(i,j,k) = s[0] * (
x(i-1,j ,k ) +
x(i+1,j ,k ))
181 + s[3] * (
x(i ,j-1,k ) +
x(i ,j+1,k ))
182 + s[5] * (
x(i ,j ,k-1) +
x(i ,j ,k+1))
183 - Real(2.)*(s[0]+s[3]+s[5]) *
x(i,j,k)
184 + Real(0.5)*s[1] * (
x(i-1,j-1,k ) +
x(i+1,j+1,k ) -
x(i-1,j+1,k ) -
x(i+1,j-1,k ))
185 + Real(0.5)*s[2] * (
x(i-1,j ,k-1) +
x(i+1,j ,k+1) -
x(i-1,j ,k+1) -
x(i+1,j ,k-1))
186 + Real(0.5)*s[4] * (
x(i ,j-1,k-1) +
x(i ,j+1,k+1) -
x(i ,j-1,k+1) -
x(i ,j+1,k-1));
198 constexpr Real omega = Real(1.25);
199 Real s0 = Real(-2.)*(s[0]+s[3]+s[5]);
200 Real Ax = s[0] * (sol(i-1,j ,k ) + sol(i+1,j ,k ))
201 + s[3] * (sol(i ,j-1,k ) + sol(i ,j+1,k ))
202 + s[5] * (sol(i ,j ,k-1) + sol(i ,j ,k+1))
204 + Real(0.5)*s[1] * (sol(i-1,j-1,k ) + sol(i+1,j+1,k ) - sol(i-1,j+1,k ) - sol(i+1,j-1,k ))
205 + Real(0.5)*s[2] * (sol(i-1,j ,k-1) + sol(i+1,j ,k+1) - sol(i-1,j ,k+1) - sol(i+1,j ,k-1))
206 + Real(0.5)*s[4] * (sol(i ,j-1,k-1) + sol(i ,j+1,k+1) - sol(i ,j-1,k+1) - sol(i ,j+1,k-1));
207 sol(i,j,k) += (rhs(i,j,k) - Ax) * (omega/s0);
211 #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
213 template <
typename HypreInt,
typename AtomicInt>
214 void mlndtslap_fill_ijmatrix_cpu (
Box const& ndbx,
215 Array4<AtomicInt const>
const& gid,
216 Array4<int const>
const& lid,
217 HypreInt*
const ncols, HypreInt*
const cols,
219 GpuArray<Real,6>
const& s) noexcept
227 HypreInt nelems_old = nelems;
229 cols[nelems] = gid(i,j,k);
230 mat[nelems] = Real(-2.)*(s[0]+s[3]+s[5]);
233 if ( gid(i ,j-1,k-1) < gidmax) {
234 cols[nelems] = gid(i ,j-1,k-1);
235 mat[nelems] = Real(0.5)*s[4];
239 if ( gid(i-1,j ,k-1) < gidmax) {
240 cols[nelems] = gid(i-1,j ,k-1);
241 mat[nelems] = Real(0.5)*s[2];
245 if ( gid(i ,j ,k-1) < gidmax) {
246 cols[nelems] = gid(i ,j ,k-1);
251 if ( gid(i+1,j ,k-1) < gidmax) {
252 cols[nelems] = gid(i+1,j ,k-1);
253 mat[nelems] = Real(-0.5)*s[2];
257 if ( gid(i ,j+1,k-1) < gidmax) {
258 cols[nelems] = gid(i ,j+1,k-1);
259 mat[nelems] = Real(-0.5)*s[4];
263 if ( gid(i-1,j-1,k ) < gidmax) {
264 cols[nelems] = gid(i-1,j-1,k );
265 mat[nelems] = Real(0.5)*s[1];
269 if ( gid(i ,j-1,k ) < gidmax) {
270 cols[nelems] = gid(i ,j-1,k );
275 if ( gid(i+1,j-1,k ) < gidmax) {
276 cols[nelems] = gid(i+1,j-1,k );
277 mat[nelems] = Real(-0.5)*s[1];
281 if ( gid(i-1,j ,k ) < gidmax) {
282 cols[nelems] = gid(i-1,j ,k );
287 if ( gid(i+1,j ,k ) < gidmax) {
288 cols[nelems] = gid(i+1,j ,k );
293 if ( gid(i-1,j+1,k ) < gidmax) {
294 cols[nelems] = gid(i-1,j+1,k );
295 mat[nelems] = Real(-0.5)*s[1];
299 if ( gid(i ,j+1,k ) < gidmax) {
300 cols[nelems] = gid(i ,j+1,k );
305 if ( gid(i+1,j+1,k ) < gidmax) {
306 cols[nelems] = gid(i+1,j+1,k );
307 mat[nelems] = Real(0.5)*s[1];
311 if ( gid(i ,j-1,k+1) < gidmax) {
312 cols[nelems] = gid(i ,j-1,k+1);
313 mat[nelems] = Real(-0.5)*s[4];
317 if ( gid(i-1,j ,k+1) < gidmax) {
318 cols[nelems] = gid(i-1,j ,k+1);
319 mat[nelems] = Real(-0.5)*s[2];
323 if ( gid(i ,j ,k+1) < gidmax) {
324 cols[nelems] = gid(i ,j ,k+1);
329 if ( gid(i+1,j ,k+1) < gidmax) {
330 cols[nelems] = gid(i+1,j ,k+1);
331 mat[nelems] = Real(0.5)*s[2];
335 if ( gid(i ,j+1,k+1) < gidmax) {
336 cols[nelems] = gid(i ,j+1,k+1);
337 mat[nelems] = Real(0.5)*s[4];
341 ncols[lid(i,j,k)] = nelems - nelems_old;
348 template <
typename HypreInt,
typename AtomicInt>
350 void mlndtslap_fill_ijmatrix_gpu (
const int ps,
const int i,
const int j,
const int k,
352 Array4<AtomicInt const>
const& gid,
353 Array4<int const>
const& lid,
354 HypreInt*
const ncols, HypreInt*
const cols,
356 GpuArray<Real,6>
const& s) noexcept
363 cols[ps] = gid(i,j,k);
364 mat[ps] = Real(-2.)*(s[0]+s[3]+s[5]);
366 if (gid(i-1,j-1,k-1) < gidmax) { ++nc; }
367 if (gid(i ,j-1,k-1) < gidmax) { ++nc; }
368 if (gid(i+1,j-1,k-1) < gidmax) { ++nc; }
369 if (gid(i-1,j ,k-1) < gidmax) { ++nc; }
370 if (gid(i ,j ,k-1) < gidmax) { ++nc; }
371 if (gid(i+1,j ,k-1) < gidmax) { ++nc; }
372 if (gid(i-1,j+1,k-1) < gidmax) { ++nc; }
373 if (gid(i ,j+1,k-1) < gidmax) { ++nc; }
374 if (gid(i+1,j+1,k-1) < gidmax) { ++nc; }
375 if (gid(i-1,j-1,k ) < gidmax) { ++nc; }
376 if (gid(i ,j-1,k ) < gidmax) { ++nc; }
377 if (gid(i+1,j-1,k ) < gidmax) { ++nc; }
378 if (gid(i-1,j ,k ) < gidmax) { ++nc; }
379 if (gid(i+1,j ,k ) < gidmax) { ++nc; }
380 if (gid(i-1,j+1,k ) < gidmax) { ++nc; }
381 if (gid(i ,j+1,k ) < gidmax) { ++nc; }
382 if (gid(i+1,j+1,k ) < gidmax) { ++nc; }
383 if (gid(i-1,j-1,k+1) < gidmax) { ++nc; }
384 if (gid(i ,j-1,k+1) < gidmax) { ++nc; }
385 if (gid(i+1,j-1,k+1) < gidmax) { ++nc; }
386 if (gid(i-1,j ,k+1) < gidmax) { ++nc; }
387 if (gid(i ,j ,k+1) < gidmax) { ++nc; }
388 if (gid(i+1,j ,k+1) < gidmax) { ++nc; }
389 if (gid(i-1,j+1,k+1) < gidmax) { ++nc; }
390 if (gid(i ,j+1,k+1) < gidmax) { ++nc; }
391 if (gid(i+1,j+1,k+1) < gidmax) { ++nc; }
392 ncols[lid(i,j,k)] = nc;
394 else if (
offset == 1 && gid(i-1,j-1,k-1) < gidmax) {
395 cols[ps] = gid(i-1,j-1,k-1);
398 else if (
offset == 2 && gid(i ,j-1,k-1) < gidmax) {
399 cols[ps] = gid(i ,j-1,k-1);
400 mat[ps] = Real(0.5)*s[4];
402 else if (
offset == 3 && gid(i+1,j-1,k-1) < gidmax) {
403 cols[ps] = gid(i+1,j-1,k-1);
406 else if (
offset == 4 && gid(i-1,j ,k-1) < gidmax) {
407 cols[ps] = gid(i-1,j ,k-1);
408 mat[ps] = Real(0.5)*s[2];
410 else if (
offset == 5 && gid(i ,j ,k-1) < gidmax) {
411 cols[ps] = gid(i ,j ,k-1);
414 else if (
offset == 6 && gid(i+1,j ,k-1) < gidmax) {
415 cols[ps] = gid(i+1,j ,k-1);
416 mat[ps] = Real(-0.5)*s[2];
418 else if (
offset == 7 && gid(i-1,j+1,k-1) < gidmax) {
419 cols[ps] = gid(i-1,j+1,k-1);
422 else if (
offset == 8 && gid(i ,j+1,k-1) < gidmax) {
423 cols[ps] = gid(i ,j+1,k-1);
424 mat[ps] = Real(-0.5)*s[4];
426 else if (
offset == 9 && gid(i+1,j+1,k-1) < gidmax) {
427 cols[ps] = gid(i+1,j+1,k-1);
430 else if (
offset == 10 && gid(i-1,j-1,k ) < gidmax) {
431 cols[ps] = gid(i-1,j-1,k );
432 mat[ps] = Real(0.5)*s[1];
434 else if (
offset == 11 && gid(i ,j-1,k ) < gidmax) {
435 cols[ps] = gid(i ,j-1,k );
438 else if (
offset == 12 && gid(i+1,j-1,k ) < gidmax) {
439 cols[ps] = gid(i+1,j-1,k );
440 mat[ps] = Real(-0.5)*s[1];
442 else if (
offset == 13 && gid(i-1,j ,k ) < gidmax) {
443 cols[ps] = gid(i-1,j ,k );
446 else if (
offset == 14 && gid(i+1,j ,k ) < gidmax) {
447 cols[ps] = gid(i+1,j ,k );
450 else if (
offset == 15 && gid(i-1,j+1,k ) < gidmax) {
451 cols[ps] = gid(i-1,j+1,k );
452 mat[ps] = Real(-0.5)*s[1];
454 else if (
offset == 16 && gid(i ,j+1,k ) < gidmax) {
455 cols[ps] = gid(i ,j+1,k );
458 else if (
offset == 17 && gid(i+1,j+1,k ) < gidmax) {
459 cols[ps] = gid(i+1,j+1,k );
460 mat[ps] = Real(0.5)*s[1];
462 else if (
offset == 18 && gid(i-1,j-1,k+1) < gidmax) {
463 cols[ps] = gid(i-1,j-1,k+1);
466 else if (
offset == 19 && gid(i ,j-1,k+1) < gidmax) {
467 cols[ps] = gid(i ,j-1,k+1);
468 mat[ps] = Real(-0.5)*s[4];
470 else if (
offset == 20 && gid(i+1,j-1,k+1) < gidmax) {
471 cols[ps] = gid(i+1,j-1,k+1);
474 else if (
offset == 21 && gid(i-1,j ,k+1) < gidmax) {
475 cols[ps] = gid(i-1,j ,k+1);
476 mat[ps] = Real(-0.5)*s[2];
478 else if (
offset == 22 && gid(i ,j ,k+1) < gidmax) {
479 cols[ps] = gid(i ,j ,k+1);
482 else if (
offset == 23 && gid(i+1,j ,k+1) < gidmax) {
483 cols[ps] = gid(i+1,j ,k+1);
484 mat[ps] = Real(0.5)*s[2];
486 else if (
offset == 24 && gid(i-1,j+1,k+1) < gidmax) {
487 cols[ps] = gid(i-1,j+1,k+1);
490 else if (
offset == 25 && gid(i ,j+1,k+1) < gidmax) {
491 cols[ps] = gid(i ,j+1,k+1);
492 mat[ps] = Real(0.5)*s[4];
494 else if (
offset == 26 && gid(i+1,j+1,k+1) < gidmax) {
495 cols[ps] = gid(i+1,j+1,k+1);
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
Array4< int const > offset
Definition: AMReX_HypreMLABecLap.cpp:1089
Array4< Real > fine
Definition: AMReX_InterpFaceRegister.cpp:90
Array4< Real const > crse
Definition: AMReX_InterpFaceRegister.cpp:92
@ max
Definition: AMReX_ParallelReduce.H:17
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ts_interp_line_x(Array4< Real const > const &crse, int ic, int jc) noexcept
Definition: AMReX_MLNodeTensorLap_2D_K.H:10
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ts_interp_line_z(Array4< Real const > const &crse, int ic, int jc, int kc) noexcept
Definition: AMReX_MLNodeTensorLap_3D_K.H:24
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ts_interp_line_y(Array4< Real const > const &crse, int ic, int jc) noexcept
Definition: AMReX_MLNodeTensorLap_2D_K.H:16
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ts_interp_face_xz(Array4< Real const > const &crse, int ic, int jc, int kc) noexcept
Definition: AMReX_MLNodeTensorLap_3D_K.H:41
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ts_interp_face_xy(Array4< Real const > const &crse, int ic, int jc) noexcept
Definition: AMReX_MLNodeTensorLap_2D_K.H:22
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ts_interp_face_yz(Array4< Real const > const &crse, int ic, int jc, int kc) noexcept
Definition: AMReX_MLNodeTensorLap_3D_K.H:51
Definition: AMReX_Amr.cpp:49
AMREX_ATTRIBUTE_FLATTEN_FOR void LoopOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:355
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndtslap_adotx(int i, int j, int k, Array4< Real > const &y, Array4< Real const > const &x, Array4< int const > const &msk, GpuArray< Real, 3 > const &s) noexcept
Definition: AMReX_MLNodeTensorLap_2D_K.H:91
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndtslap_gauss_seidel(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< int const > const &msk, GpuArray< Real, 3 > const &s) noexcept
Definition: AMReX_MLNodeTensorLap_2D_K.H:105
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndtslap_semi_interpadd(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, int) noexcept
Definition: AMReX_MLNodeTensorLap_1D_K.H:13
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndtslap_interpadd(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition: AMReX_MLNodeTensorLap_1D_K.H:8
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo/rati...
Definition: AMReX_Box.H:1304
Definition: AMReX_Array4.H:61
Definition: AMReX_Array.H:34