Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLALap_3D_K.H
Go to the documentation of this file.
1#ifndef AMREX_MLALAP_3D_K_H_
2#define AMREX_MLALAP_3D_K_H_
3#include <AMReX_Config.H>
4
5namespace amrex {
6
7template <typename RT>
9void mlalap_adotx (Box const& box, Array4<RT> const& y,
10 Array4<RT const> const& x,
11 Array4<RT const> const& a,
12 GpuArray<RT,AMREX_SPACEDIM> const& dxinv,
13 RT alpha, RT beta, int ncomp) noexcept
14{
15 const RT dhx = beta*dxinv[0]*dxinv[0];
16 const RT dhy = beta*dxinv[1]*dxinv[1];
17 const RT dhz = beta*dxinv[2]*dxinv[2];
18
19 const auto lo = amrex::lbound(box);
20 const auto hi = amrex::ubound(box);
21
22 for (int n = 0; n < ncomp; ++n) {
23 for (int k = lo.z; k <= hi.z; ++k) {
24 for (int j = lo.y; j <= hi.y; ++j) {
26 for (int i = lo.x; i <= hi.x; ++i) {
27 y(i,j,k,n) = alpha*a(i,j,k,n)*x(i,j,k,n)
28 - dhx * (x(i-1,j,k,n) - RT(2.0)*x(i,j,k,n) + x(i+1,j,k,n))
29 - dhy * (x(i,j-1,k,n) - RT(2.0)*x(i,j,k,n) + x(i,j+1,k,n))
30 - dhz * (x(i,j,k-1,n) - RT(2.0)*x(i,j,k,n) + x(i,j,k+1,n));
31 }
32 }
33 }
34 }
35}
36
37template <typename RT>
39void mlalap_normalize (Box const& box, Array4<RT> const& x,
40 Array4<RT const> const& a,
41 GpuArray<RT,AMREX_SPACEDIM> const& dxinv,
42 RT alpha, RT beta, int ncomp) noexcept
43{
44 const RT dhx = beta*dxinv[0]*dxinv[0];
45 const RT dhy = beta*dxinv[1]*dxinv[1];
46 const RT dhz = beta*dxinv[2]*dxinv[2];
47 const RT fac = RT(2.0)*(dhx+dhy+dhz);
48
49 const auto lo = amrex::lbound(box);
50 const auto hi = amrex::ubound(box);
51
52 for (int n = 0; n < ncomp; ++n) {
53 for (int k = lo.z; k <= hi.z; ++k) {
54 for (int j = lo.y; j <= hi.y; ++j) {
56 for (int i = lo.x; i <= hi.x; ++i) {
57 x(i,j,k,n) /= alpha*a(i,j,k,n) + fac;
58 }
59 }
60 }
61 }
62}
63
64template <typename RT>
66void mlalap_flux_x (Box const& box, Array4<RT> const& fx, Array4<RT const> const& sol,
67 RT fac, int ncomp) noexcept
68{
69 const auto lo = amrex::lbound(box);
70 const auto hi = amrex::ubound(box);
71
72 for (int n = 0; n < ncomp; ++n) {
73 for (int k = lo.z; k <= hi.z; ++k) {
74 for (int j = lo.y; j <= hi.y; ++j) {
76 for (int i = lo.x; i <= hi.x; ++i) {
77 fx(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i-1,j,k,n));
78 }
79 }
80 }
81 }
82}
83
84template <typename RT>
86void mlalap_flux_xface (Box const& box, Array4<RT> const& fx, Array4<RT const> const& sol,
87 RT fac, int xlen, int ncomp) noexcept
88{
89 const auto lo = amrex::lbound(box);
90 const auto hi = amrex::ubound(box);
91
92 for (int n = 0; n < ncomp; ++n) {
93 for (int k = lo.z; k <= hi.z; ++k) {
94 for (int j = lo.y; j <= hi.y; ++j) {
95 int i = lo.x;
96 fx(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i-1,j,k,n));
97 i += xlen;
98 fx(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i-1,j,k,n));
99 }
100 }
101 }
102}
103
104template <typename RT>
106void mlalap_flux_y (Box const& box, Array4<RT> const& fy, Array4<RT const> const& sol,
107 RT fac, int ncomp) noexcept
108{
109 const auto lo = amrex::lbound(box);
110 const auto hi = amrex::ubound(box);
111
112 for (int n = 0; n < ncomp; ++n) {
113 for (int k = lo.z; k <= hi.z; ++k) {
114 for (int j = lo.y; j <= hi.y; ++j) {
116 for (int i = lo.x; i <= hi.x; ++i) {
117 fy(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i,j-1,k,n));
118 }
119 }
120 }
121 }
122}
123
124template <typename RT>
126void mlalap_flux_yface (Box const& box, Array4<RT> const& fy, Array4<RT const> const& sol,
127 RT fac, int ylen, int ncomp) noexcept
128{
129 const auto lo = amrex::lbound(box);
130 const auto hi = amrex::ubound(box);
131
132 for (int n = 0; n < ncomp; ++n) {
133 for (int k = lo.z; k <= hi.z; ++k) {
134 int j = lo.y;
136 for (int i = lo.x; i <= hi.x; ++i) {
137 fy(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i,j-1,k,n));
138 }
139 j += ylen;
141 for (int i = lo.x; i <= hi.x; ++i) {
142 fy(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i,j-1,k,n));
143 }
144 }
145 }
146}
147
148template <typename RT>
150void mlalap_flux_z (Box const& box, Array4<RT> const& fz, Array4<RT const> const& sol,
151 RT fac, int ncomp) noexcept
152{
153 const auto lo = amrex::lbound(box);
154 const auto hi = amrex::ubound(box);
155
156 for (int n = 0; n < ncomp; ++n) {
157 for (int k = lo.z; k <= hi.z; ++k) {
158 for (int j = lo.y; j <= hi.y; ++j) {
160 for (int i = lo.x; i <= hi.x; ++i) {
161 fz(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i,j,k-1,n));
162 }
163 }
164 }
165 }
166}
167
168template <typename RT>
170void mlalap_flux_zface (Box const& box, Array4<RT> const& fz, Array4<RT const> const& sol,
171 RT fac, int zlen, int ncomp) noexcept
172{
173 const auto lo = amrex::lbound(box);
174 const auto hi = amrex::ubound(box);
175
176 for (int n = 0; n < ncomp; ++n) {
177 int k = lo.z;
178 for (int j = lo.y; j <= hi.y; ++j) {
180 for (int i = lo.x; i <= hi.x; ++i) {
181 fz(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i,j,k-1,n));
182 }
183 }
184
185 k += zlen;
186 for (int j = lo.y; j <= hi.y; ++j) {
188 for (int i = lo.x; i <= hi.x; ++i) {
189 fz(i,j,k,n) = -fac*(sol(i,j,k,n)-sol(i,j,k-1,n));
190 }
191 }
192 }
193}
194
195template <typename RT>
197void mlalap_gsrb (Box const& box, Array4<RT> const& phi,
198 Array4<RT const> const& rhs, RT alpha,
199 RT dhx, RT dhy, RT dhz, Array4<RT const> const& a,
200 Array4<RT const> const& f0, Array4<int const> const& m0,
201 Array4<RT const> const& f1, Array4<int const> const& m1,
202 Array4<RT const> const& f2, Array4<int const> const& m2,
203 Array4<RT const> const& f3, Array4<int const> const& m3,
204 Array4<RT const> const& f4, Array4<int const> const& m4,
205 Array4<RT const> const& f5, Array4<int const> const& m5,
206 Box const& vbox, int redblack, int ncomp) noexcept
207{
208 const auto lo = amrex::lbound(box);
209 const auto hi = amrex::ubound(box);
210 const auto vlo = amrex::lbound(vbox);
211 const auto vhi = amrex::ubound(vbox);
212
213 constexpr RT omega = RT(1.15);
214
215 const RT dhfac = RT(2.)*(dhx+dhy+dhz);
216
217 for (int n = 0; n < ncomp; ++n) {
218 for (int k = lo.z; k <= hi.z; ++k) {
219 for (int j = lo.y; j <= hi.y; ++j) {
221 for (int i = lo.x; i <= hi.x; ++i) {
222 if ((i+j+k+redblack)%2 == 0) {
223 RT cf0 = (i == vlo.x && m0(vlo.x-1,j,k) > 0)
224 ? f0(vlo.x,j,k,n) : RT(0.0);
225 RT cf1 = (j == vlo.y && m1(i,vlo.y-1,k) > 0)
226 ? f1(i,vlo.y,k,n) : RT(0.0);
227 RT cf2 = (k == vlo.z && m2(i,j,vlo.z-1) > 0)
228 ? f2(i,j,vlo.z,n) : RT(0.0);
229 RT cf3 = (i == vhi.x && m3(vhi.x+1,j,k) > 0)
230 ? f3(vhi.x,j,k,n) : RT(0.0);
231 RT cf4 = (j == vhi.y && m4(i,vhi.y+1,k) > 0)
232 ? f4(i,vhi.y,k,n) : RT(0.0);
233 RT cf5 = (k == vhi.z && m5(i,j,vhi.z+1) > 0)
234 ? f5(i,j,vhi.z,n) : RT(0.0);
235
236 RT gamma = alpha*a(i,j,k,n) + dhfac;
237
238 RT g_m_d = gamma - dhx*(cf0+cf3) - dhy*(cf1+cf4) - dhz*(cf2+cf5);
239
240 RT rho = dhx*(phi(i-1,j,k,n) + phi(i+1,j,k,n))
241 + dhy*(phi(i,j-1,k,n) + phi(i,j+1,k,n))
242 + dhz*(phi(i,j,k-1,n) + phi(i,j,k+1,n));
243
244 RT res = rhs(i,j,k,n) - (gamma*phi(i,j,k,n) - rho);
245 phi(i,j,k,n) = phi(i,j,k,n) + omega/g_m_d * res;
246 }
247 }
248 }
249 }
250 }
251}
252
253}
254#endif
#define AMREX_PRAGMA_SIMD
Definition AMReX_Extension.H:80
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_x(Box const &box, Array4< RT > const &fx, Array4< RT const > const &sol, RT fac, int ncomp) noexcept
Definition AMReX_MLALap_1D_K.H:101
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_normalize(Box const &box, Array4< RT > const &x, Array4< RT const > const &a, GpuArray< RT, AMREX_SPACEDIM > const &dxinv, RT alpha, RT beta, int ncomp) noexcept
Definition AMReX_MLALap_1D_K.H:58
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_yface(Box const &box, Array4< RT > const &fy, Array4< RT const > const &sol, RT fac, int ylen, int ncomp) noexcept
Definition AMReX_MLALap_3D_K.H:126
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_gsrb(Box const &box, Array4< RT > const &phi, Array4< RT const > const &rhs, RT alpha, RT dhx, Array4< RT const > const &a, Array4< RT const > const &f0, Array4< int const > const &m0, Array4< RT const > const &f1, Array4< int const > const &m1, Box const &vbox, int redblack, int ncomp) noexcept
Definition AMReX_MLALap_1D_K.H:168
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_adotx(Box const &box, Array4< RT > const &y, Array4< RT const > const &x, Array4< RT const > const &a, GpuArray< RT, AMREX_SPACEDIM > const &dxinv, RT alpha, RT beta, int ncomp) noexcept
Definition AMReX_MLALap_1D_K.H:9
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 Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:308
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_zface(Box const &box, Array4< RT > const &fz, Array4< RT const > const &sol, RT fac, int zlen, int ncomp) noexcept
Definition AMReX_MLALap_3D_K.H:170
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_y(Box const &box, Array4< RT > const &fy, Array4< RT const > const &sol, RT fac, int ncomp) noexcept
Definition AMReX_MLALap_3D_K.H:106
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_xface(Box const &box, Array4< RT > const &fx, Array4< RT const > const &sol, RT fac, int xlen, int ncomp) noexcept
Definition AMReX_MLALap_1D_K.H:135
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_z(Box const &box, Array4< RT > const &fz, Array4< RT const > const &sol, RT fac, int ncomp) noexcept
Definition AMReX_MLALap_3D_K.H:150
Definition AMReX_Array4.H:61