Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLALap_2D_K.H
Go to the documentation of this file.
1#ifndef AMREX_MLALAP_2D_K_H_
2#define AMREX_MLALAP_2D_K_H_
3#include <AMReX_Config.H>
4
5namespace amrex::TwoD {
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,2> 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
18 const auto lo = amrex::lbound(box);
19 const auto hi = amrex::ubound(box);
20
21 for (int n = 0; n < ncomp; ++n) {
22 for (int j = lo.y; j <= hi.y; ++j) {
24 for (int i = lo.x; i <= hi.x; ++i) {
25 y(i,j,0,n) = alpha*a(i,j,0,n)*x(i,j,0,n)
26 - dhx * (x(i-1,j,0,n) - RT(2.)*x(i,j,0,n) + x(i+1,j,0,n))
27 - dhy * (x(i,j-1,0,n) - RT(2.)*x(i,j,0,n) + x(i,j+1,0,n));
28 }
29 }
30 }
31}
32
33template <typename RT>
35void mlalap_adotx_m (Box const& box, Array4<RT> const& y,
36 Array4<RT const> const& x,
37 Array4<RT const> const& a,
38 GpuArray<RT,2> const& dxinv,
39 RT alpha, RT beta, RT dx, RT probxlo, int ncomp) noexcept
40{
41 const RT dhx = beta*dxinv[0]*dxinv[0];
42 const RT dhy = beta*dxinv[1]*dxinv[1];
43
44 const auto lo = amrex::lbound(box);
45 const auto hi = amrex::ubound(box);
46
47 for (int n = 0; n < ncomp; ++n) {
48 for (int j = lo.y; j <= hi.y; ++j) {
50 for (int i = lo.x; i <= hi.x; ++i) {
51 RT rel = probxlo + i*dx;
52 RT rer = probxlo +(i+1)*dx;
53 RT rc = probxlo + (i+RT(0.5))*dx;
54 y(i,j,0,n) = alpha*a(i,j,0,n)*x(i,j,0,n)*rc
55 - dhx * (rer*(x(i+1,j,0,n) - x(i ,j,0,n))
56 - rel*(x(i ,j,0,n) - x(i-1,j,0,n)))
57 - dhy * rc * (x(i,j-1,0,n) - RT(2.)*x(i,j,0,n) + x(i,j+1,0,n));
58 }
59 }
60 }
61}
62
63template <typename RT>
65void mlalap_normalize (Box const& box, Array4<RT> const& x,
66 Array4<RT const> const& a,
67 GpuArray<RT,2> const& dxinv,
68 RT alpha, RT beta, int ncomp) noexcept
69{
70 const RT dhx = beta*dxinv[0]*dxinv[0];
71 const RT dhy = beta*dxinv[1]*dxinv[1];
72
73 const auto lo = amrex::lbound(box);
74 const auto hi = amrex::ubound(box);
75
76 for (int n = 0; n < ncomp; ++n) {
77 for (int j = lo.y; j <= hi.y; ++j) {
79 for (int i = lo.x; i <= hi.x; ++i) {
80 x(i,j,0,n) /= alpha*a(i,j,0,n) + RT(2.0)*(dhx + dhy);
81 }
82 }
83 }
84}
85
86template <typename RT>
88void mlalap_normalize_m (Box const& box, Array4<RT> const& x,
89 Array4<RT const> const& a,
90 GpuArray<RT,2> const& dxinv,
91 RT alpha, RT beta,
92 RT dx, RT probxlo, int ncomp) noexcept
93{
94 const RT dhx = beta*dxinv[0]*dxinv[0];
95 const RT dhy = beta*dxinv[1]*dxinv[1];
96
97 const auto lo = amrex::lbound(box);
98 const auto hi = amrex::ubound(box);
99
100 for (int n = 0; n < ncomp; ++n) {
101 for (int j = lo.y; j <= hi.y; ++j) {
103 for (int i = lo.x; i <= hi.x; ++i) {
104 RT rel = probxlo + i*dx;
105 RT rer = probxlo +(i+1)*dx;
106 RT rc = probxlo + (i+RT(0.5))*dx;
107 x(i,j,0,n) /= alpha*a(i,j,0,n)*rc
108 + dhx*(rel+rer) + dhy*(rc*RT(2.0));
109 }
110 }
111 }
112}
113
114template <typename RT>
116void mlalap_flux_x (Box const& box, Array4<RT> const& fx,
117 Array4<RT const> const& sol, RT fac, int ncomp) noexcept
118{
119 const auto lo = amrex::lbound(box);
120 const auto hi = amrex::ubound(box);
121
122 for (int n = 0; n < ncomp; ++n) {
123 for (int j = lo.y; j <= hi.y; ++j) {
125 for (int i = lo.x; i <= hi.x; ++i) {
126 fx(i,j,0,n) = -fac*(sol(i,j,0,n)-sol(i-1,j,0,n));
127 }
128 }
129 }
130}
131
132template <typename RT>
134void mlalap_flux_x_m (Box const& box, Array4<RT> const& fx,
135 Array4<RT const> const& sol, RT fac,
136 RT dx, RT probxlo, int ncomp) noexcept
137{
138 const auto lo = amrex::lbound(box);
139 const auto hi = amrex::ubound(box);
140
141 for (int n = 0; n < ncomp; ++n) {
142 for (int j = lo.y; j <= hi.y; ++j) {
144 for (int i = lo.x; i <= hi.x; ++i) {
145 RT re = probxlo + i*dx;
146 fx(i,j,0,n) = -fac*re*(sol(i,j,0,n)-sol(i-1,j,0,n));
147 }
148 }
149 }
150}
151
152template <typename RT>
154void mlalap_flux_xface (Box const& box, Array4<RT> const& fx,
155 Array4<RT const> const& sol, RT fac, int xlen, int ncomp) noexcept
156{
157 const auto lo = amrex::lbound(box);
158 const auto hi = amrex::ubound(box);
159
160 for (int n = 0; n < ncomp; ++n) {
161 for (int j = lo.y; j <= hi.y; ++j) {
162 int i = lo.x;
163 fx(i,j,0,n) = -fac*(sol(i,j,0,n)-sol(i-1,j,0,n));
164 i += xlen;
165 fx(i,j,0,n) = -fac*(sol(i,j,0,n)-sol(i-1,j,0,n));
166 }
167 }
168}
169
170template <typename RT>
172void mlalap_flux_xface_m (Box const& box, Array4<RT> const& fx,
173 Array4<RT const> const& sol, RT fac, int xlen,
174 RT dx, RT probxlo, int ncomp) noexcept
175{
176 const auto lo = amrex::lbound(box);
177 const auto hi = amrex::ubound(box);
178
179 for (int n = 0; n < ncomp; ++n) {
180 for (int j = lo.y; j <= hi.y; ++j) {
181 int i = lo.x;
182 RT re = probxlo + i*dx;
183 fx(i,j,0,n) = -fac*re*(sol(i,j,0,n)-sol(i-1,j,0,n));
184 i += xlen;
185 re = probxlo + i*dx;
186 fx(i,j,0,n) = -fac*re*(sol(i,j,0,n)-sol(i-1,j,0,n));
187 }
188 }
189}
190
191template <typename RT>
193void mlalap_flux_y (Box const& box, Array4<RT> const& fy,
194 Array4<RT const> const& sol, RT fac, int ncomp) noexcept
195{
196 const auto lo = amrex::lbound(box);
197 const auto hi = amrex::ubound(box);
198
199 for (int n = 0; n < ncomp; ++n) {
200 for (int j = lo.y; j <= hi.y; ++j) {
202 for (int i = lo.x; i <= hi.x; ++i) {
203 fy(i,j,0,n) = -fac*(sol(i,j,0,n)-sol(i,j-1,0,n));
204 }
205 }
206 }
207}
208
209template <typename RT>
211void mlalap_flux_y_m (Box const& box, Array4<RT> const& fy,
212 Array4<RT const> const& sol, RT fac,
213 RT dx, RT probxlo, int ncomp) noexcept
214{
215 const auto lo = amrex::lbound(box);
216 const auto hi = amrex::ubound(box);
217
218 for (int n = 0; n < ncomp; ++n) {
219 for (int j = lo.y; j <= hi.y; ++j) {
221 for (int i = lo.x; i <= hi.x; ++i) {
222 RT rc = probxlo + (i+RT(0.5))*dx;
223 fy(i,j,0,n) = -fac*rc*(sol(i,j,0,n)-sol(i,j-1,0,n));
224 }
225 }
226 }
227}
228
229template <typename RT>
231void mlalap_flux_yface (Box const& box, Array4<RT> const& fy,
232 Array4<RT const> const& sol, RT fac, int ylen, int ncomp) noexcept
233{
234 const auto lo = amrex::lbound(box);
235 const auto hi = amrex::ubound(box);
236
237 for (int n = 0; n < ncomp; ++n) {
238 int j = lo.y;
240 for (int i = lo.x; i <= hi.x; ++i) {
241 fy(i,j,0,n) = -fac*(sol(i,j,0,n)-sol(i,j-1,0,n));
242 }
243 j += ylen;
245 for (int i = lo.x; i <= hi.x; ++i) {
246 fy(i,j,0,n) = -fac*(sol(i,j,0,n)-sol(i,j-1,0,n));
247 }
248 }
249}
250
251template <typename RT>
253void mlalap_flux_yface_m (Box const& box, Array4<RT> const& fy,
254 Array4<RT const> const& sol, RT fac, int ylen,
255 RT dx, RT probxlo, int ncomp) noexcept
256{
257 const auto lo = amrex::lbound(box);
258 const auto hi = amrex::ubound(box);
259
260 for (int n = 0; n < ncomp; ++n) {
261 int j = lo.y;
263 for (int i = lo.x; i <= hi.x; ++i) {
264 RT rc = probxlo + (i+RT(0.5))*dx;
265 fy(i,j,0,n) = -fac*rc*(sol(i,j,0,n)-sol(i,j-1,0,n));
266 }
267 j += ylen;
269 for (int i = lo.x; i <= hi.x; ++i) {
270 RT rc = probxlo + (i+RT(0.5))*dx;
271 fy(i,j,0,n) = -fac*rc*(sol(i,j,0,n)-sol(i,j-1,0,n));
272 }
273 }
274}
275
276template <typename RT>
278void mlalap_gsrb (Box const& box, Array4<RT> const& phi,
279 Array4<RT const> const& rhs, RT alpha,
280 RT dhx, RT dhy, Array4<RT const> const& a,
281 Array4<RT const> const& f0, Array4<int const> const& m0,
282 Array4<RT const> const& f1, Array4<int const> const& m1,
283 Array4<RT const> const& f2, Array4<int const> const& m2,
284 Array4<RT const> const& f3, Array4<int const> const& m3,
285 Box const& vbox, int redblack, int ncomp) noexcept
286{
287 const auto lo = amrex::lbound(box);
288 const auto hi = amrex::ubound(box);
289 const auto vlo = amrex::lbound(vbox);
290 const auto vhi = amrex::ubound(vbox);
291
292 for (int n = 0; n < ncomp; ++n) {
293 for (int j = lo.y; j <= hi.y; ++j) {
295 for (int i = lo.x; i <= hi.x; ++i) {
296 if ((i+j+redblack)%2 == 0) {
297 RT cf0 = (i == vlo.x && m0(vlo.x-1,j,0) > 0)
298 ? f0(vlo.x,j,0,n) : RT(0.0);
299 RT cf1 = (j == vlo.y && m1(i,vlo.y-1,0) > 0)
300 ? f1(i,vlo.y,0,n) : RT(0.0);
301 RT cf2 = (i == vhi.x && m2(vhi.x+1,j,0) > 0)
302 ? f2(vhi.x,j,0,n) : RT(0.0);
303 RT cf3 = (j == vhi.y && m3(i,vhi.y+1,0) > 0)
304 ? f3(i,vhi.y,0,n) : RT(0.0);
305
306 RT delta = dhx*(cf0 + cf2) + dhy*(cf1 + cf3);
307
308 RT gamma = alpha*a(i,j,0,n) + RT(2.0)*(dhx + dhy);
309
310 RT rho = dhx*(phi(i-1,j,0,n) + phi(i+1,j,0,n))
311 + dhy*(phi(i,j-1,0,n) + phi(i,j+1,0,n));
312
313 phi(i,j,0,n) = (rhs(i,j,0,n) + rho - phi(i,j,0,n)*delta)
314 / (gamma - delta);
315 }
316 }
317 }
318 }
319}
320
321template <typename RT>
323void mlalap_gsrb_m (Box const& box, Array4<RT> const& phi,
324 Array4<RT const> const& rhs, RT alpha,
325 RT dhx, RT dhy, Array4<RT const> const& a,
326 Array4<RT const> const& f0, Array4<int const> const& m0,
327 Array4<RT const> const& f1, Array4<int const> const& m1,
328 Array4<RT const> const& f2, Array4<int const> const& m2,
329 Array4<RT const> const& f3, Array4<int const> const& m3,
330 Box const& vbox, int redblack,
331 RT dx, RT probxlo, int ncomp) noexcept
332{
333 const auto lo = amrex::lbound(box);
334 const auto hi = amrex::ubound(box);
335 const auto vlo = amrex::lbound(vbox);
336 const auto vhi = amrex::ubound(vbox);
337
338 for (int n = 0; n < ncomp; ++n) {
339 for (int j = lo.y; j <= hi.y; ++j) {
341 for (int i = lo.x; i <= hi.x; ++i) {
342 if ((i+j+redblack)%2 == 0) {
343 RT cf0 = (i == vlo.x && m0(vlo.x-1,j,0) > 0)
344 ? f0(vlo.x,j,0,n) : RT(0.0);
345 RT cf1 = (j == vlo.y && m1(i,vlo.y-1,0) > 0)
346 ? f1(i,vlo.y,0,n) : RT(0.0);
347 RT cf2 = (i == vhi.x && m2(vhi.x+1,j,0) > 0)
348 ? f2(vhi.x,j,0,n) : RT(0.0);
349 RT cf3 = (j == vhi.y && m3(i,vhi.y+1,0) > 0)
350 ? f3(i,vhi.y,0,n) : RT(0.0);
351
352 RT rel = probxlo + i*dx;
353 RT rer = probxlo +(i+1)*dx;
354 RT rc = probxlo + (i+RT(0.5))*dx;
355
356 RT delta = dhx*(rel*cf0 + rer*cf2) + dhy*rc*(cf1 + cf3);
357
358 RT gamma = alpha*a(i,j,0,n)*rc + dhx*(rel+rer) + dhy*(RT(2.)*rc);
359
360 RT rho = dhx*(rel*phi(i-1,j,0,n) + rer*phi(i+1,j,0,n))
361 + dhy* rc*(phi(i,j-1,0,n) + phi(i,j+1,0,n));
362
363 phi(i,j,0,n) = (rhs(i,j,0,n) + rho - phi(i,j,0,n)*delta)
364 / (gamma - delta);
365 }
366 }
367 }
368 }
369}
370
371}
372
373#if (AMREX_SPACEDIM == 2)
374namespace amrex {
375 using namespace TwoD;
376}
377#endif
378
379#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_MLALap_2D_K.H:5
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_2D_K.H:116
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_normalize_m(Box const &box, Array4< RT > const &x, Array4< RT const > const &a, GpuArray< RT, 2 > const &dxinv, RT alpha, RT beta, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:88
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_2D_K.H:154
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_2D_K.H:193
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, 2 > const &dxinv, RT alpha, RT beta, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:9
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_xface_m(Box const &box, Array4< RT > const &fx, Array4< RT const > const &sol, RT fac, int xlen, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:172
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_x_m(Box const &box, Array4< RT > const &fx, Array4< RT const > const &sol, RT fac, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:134
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, RT dhy, 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, Array4< RT const > const &f2, Array4< int const > const &m2, Array4< RT const > const &f3, Array4< int const > const &m3, Box const &vbox, int redblack, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:278
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_yface_m(Box const &box, Array4< RT > const &fy, Array4< RT const > const &sol, RT fac, int ylen, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:253
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_gsrb_m(Box const &box, Array4< RT > const &phi, Array4< RT const > const &rhs, RT alpha, RT dhx, RT dhy, 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, Array4< RT const > const &f2, Array4< int const > const &m2, Array4< RT const > const &f3, Array4< int const > const &m3, Box const &vbox, int redblack, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:323
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_flux_y_m(Box const &box, Array4< RT > const &fy, Array4< RT const > const &sol, RT fac, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:211
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_adotx_m(Box const &box, Array4< RT > const &y, Array4< RT const > const &x, Array4< RT const > const &a, GpuArray< RT, 2 > const &dxinv, RT alpha, RT beta, RT dx, RT probxlo, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:35
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlalap_normalize(Box const &box, Array4< RT > const &x, Array4< RT const > const &a, GpuArray< RT, 2 > const &dxinv, RT alpha, RT beta, int ncomp) noexcept
Definition AMReX_MLALap_2D_K.H:65
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_2D_K.H:231
Definition AMReX_Amr.cpp:49
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
Definition AMReX_Array4.H:61
Definition AMReX_Array.H:34