Block-Structured AMR Software Framework
AMReX_MLNodeABecLap_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_MLNODEABECLAP_3D_K_H_
2 #define AMREX_MLNODEABECLAP_3D_K_H_
3 
4 namespace amrex {
5 
6 inline void
7 mlndabeclap_gauss_seidel_aa (Box const& bx, Array4<Real> const& sol,
8  Array4<Real const> const& rhs,
9  Real alpha, Real beta,
10  Array4<Real const> const& acf,
11  Array4<Real const> const& bcf,
12  Array4<int const> const& msk,
13  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
14 {
15  Real facx = Real(1.0/36.0)*dxinv[0]*dxinv[0];
16  Real facy = Real(1.0/36.0)*dxinv[1]*dxinv[1];
17  Real facz = Real(1.0/36.0)*dxinv[2]*dxinv[2];
18  Real fxyz = facx + facy + facz;
19  Real fmx2y2z = -facx + Real(2.0)*facy + Real(2.0)*facz;
20  Real f2xmy2z = Real(2.0)*facx - facy + Real(2.0)*facz;
21  Real f2x2ymz = Real(2.0)*facx + Real(2.0)*facy - facz;
22  Real f4xm2ym2z = Real(4.0)*facx - Real(2.0)*facy - Real(2.0)*facz;
23  Real fm2x4ym2z = -Real(2.0)*facx + Real(4.0)*facy - Real(2.0)*facz;
24  Real fm2xm2y4z = -Real(2.0)*facx - Real(2.0)*facy + Real(4.0)*facz;
25 
26  amrex::LoopOnCpu(bx, [=] (int i, int j, int k) noexcept
27  {
28  if (msk(i,j,k)) {
29  sol(i,j,k) = Real(0.0);
30  } else {
31  Real s0 = Real(-4.0)*fxyz*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i,j,k-1)
32  +bcf(i-1,j-1,k )+bcf(i,j-1,k )+bcf(i-1,j,k )+bcf(i,j,k ));
33  Real lap = sol(i,j,k)*s0
34  + fxyz*(sol(i-1,j-1,k-1)*bcf(i-1,j-1,k-1)
35  + sol(i+1,j-1,k-1)*bcf(i ,j-1,k-1)
36  + sol(i-1,j+1,k-1)*bcf(i-1,j ,k-1)
37  + sol(i+1,j+1,k-1)*bcf(i ,j ,k-1)
38  + sol(i-1,j-1,k+1)*bcf(i-1,j-1,k )
39  + sol(i+1,j-1,k+1)*bcf(i ,j-1,k )
40  + sol(i-1,j+1,k+1)*bcf(i-1,j ,k )
41  + sol(i+1,j+1,k+1)*bcf(i ,j ,k ))
42  + fmx2y2z*(sol(i ,j-1,k-1)*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1))
43  + sol(i ,j+1,k-1)*(bcf(i-1,j ,k-1)+bcf(i,j ,k-1))
44  + sol(i ,j-1,k+1)*(bcf(i-1,j-1,k )+bcf(i,j-1,k ))
45  + sol(i ,j+1,k+1)*(bcf(i-1,j ,k )+bcf(i,j ,k )))
46  + f2xmy2z*(sol(i-1,j ,k-1)*(bcf(i-1,j-1,k-1)+bcf(i-1,j,k-1))
47  + sol(i+1,j ,k-1)*(bcf(i ,j-1,k-1)+bcf(i ,j,k-1))
48  + sol(i-1,j ,k+1)*(bcf(i-1,j-1,k )+bcf(i-1,j,k ))
49  + sol(i+1,j ,k+1)*(bcf(i ,j-1,k )+bcf(i ,j,k )))
50  + f2x2ymz*(sol(i-1,j-1,k )*(bcf(i-1,j-1,k-1)+bcf(i-1,j-1,k))
51  + sol(i+1,j-1,k )*(bcf(i ,j-1,k-1)+bcf(i ,j-1,k))
52  + sol(i-1,j+1,k )*(bcf(i-1,j ,k-1)+bcf(i-1,j ,k))
53  + sol(i+1,j+1,k )*(bcf(i ,j ,k-1)+bcf(i ,j ,k)))
54  + f4xm2ym2z*(sol(i-1,j,k)*(bcf(i-1,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i-1,j-1,k)+bcf(i-1,j,k))
55  + sol(i+1,j,k)*(bcf(i ,j-1,k-1)+bcf(i ,j,k-1)+bcf(i ,j-1,k)+bcf(i ,j,k)))
56  + fm2x4ym2z*(sol(i,j-1,k)*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j-1,k)+bcf(i,j-1,k))
57  + sol(i,j+1,k)*(bcf(i-1,j ,k-1)+bcf(i,j ,k-1)+bcf(i-1,j ,k)+bcf(i,j ,k)))
58  + fm2xm2y4z*(sol(i,j,k-1)*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i,j,k-1))
59  + sol(i,j,k+1)*(bcf(i-1,j-1,k )+bcf(i,j-1,k )+bcf(i-1,j,k )+bcf(i,j,k )));
60  Real Ax = alpha*acf(i,j,k)*sol(i,j,k) - beta*lap;
61 
62  sol(i,j,k) += (rhs(i,j,k) - Ax) / (alpha*acf(i,j,k)-beta*s0);
63  }
64  });
65 }
66 
67 AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
68 mlndabeclap_jacobi_aa (int i, int j, int k, Array4<Real> const& sol,
69  Real lap, Array4<Real const> const& rhs,
70  Real alpha, Real beta,
71  Array4<Real const> const& acf,
72  Array4<Real const> const& bcf,
73  Array4<int const> const& msk,
74  GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
75 {
76  if (msk(i,j,k)) {
77  sol(i,j,k) = Real(0.0);
78  } else {
79  Real fxyz = Real(-4.0 / 36.0)*(dxinv[0]*dxinv[0] +
80  dxinv[1]*dxinv[1] +
81  dxinv[2]*dxinv[2]);
82  Real s0 = fxyz*(bcf(i-1,j-1,k-1)+bcf(i,j-1,k-1)+bcf(i-1,j,k-1)+bcf(i,j,k-1)
83  +bcf(i-1,j-1,k )+bcf(i,j-1,k )+bcf(i-1,j,k )+bcf(i,j,k));
84  Real Ax = alpha*acf(i,j,k)*sol(i,j,k) - beta*lap;
85 
86  sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
87  / (alpha*acf(i,j,k)-beta*s0);
88  }
89 }
90 
91 }
92 
93 #endif
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
amrex
Definition: AMReX_Amr.cpp:49
amrex::LoopOnCpu
AMREX_ATTRIBUTE_FLATTEN_FOR void LoopOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition: AMReX_Loop.H:355
amrex::Box
BoxND< AMREX_SPACEDIM > Box
Definition: AMReX_BaseFwd.H:27
amrex::mlndabeclap_jacobi_aa
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndabeclap_jacobi_aa(int, int, int, Array4< Real > const &, Real, Array4< Real const > const &, Real, Real, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, GpuArray< Real, AMREX_SPACEDIM > const &) noexcept
Definition: AMReX_MLNodeABecLap_1D_K.H:17
amrex::mlndabeclap_gauss_seidel_aa
void mlndabeclap_gauss_seidel_aa(Box const &, Array4< Real > const &, Array4< Real const > const &, Real, Real, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, GpuArray< Real, AMREX_SPACEDIM > const &) noexcept
Definition: AMReX_MLNodeABecLap_1D_K.H:7