Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLNodeLap_2D_K.H
Go to the documentation of this file.
1#ifndef AMREX_MLNODELAP_2D_K_H_
2#define AMREX_MLNODELAP_2D_K_H_
3#include <AMReX_Config.H>
4
5namespace amrex {
6
7AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
8void mlndlap_zero_fine (int i, int j, int, Array4<Real> const& phi,
9 Array4<int const> const& msk, int fine_flag) noexcept
10{
11 // Testing if the node is covered by a fine level in computing
12 // coarse sync residual
13 if (msk(i-1,j-1,0) == fine_flag &&
14 msk(i ,j-1,0) == fine_flag &&
15 msk(i-1,j ,0) == fine_flag &&
16 msk(i ,j ,0) == fine_flag)
17 {
18 phi(i,j,0) = Real(0.0);
19 }
20}
21
22//
23// coeffs
24//
25
26AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
27void mlndlap_avgdown_coeff_x (int i, int j, int k, Array4<Real> const& crse,
28 Array4<Real const> const& fine) noexcept
29{
30 Real a = fine(2*i ,2*j,k) + fine(2*i ,2*j+1,k);
31 Real b = fine(2*i+1,2*j,k) + fine(2*i+1,2*j+1,k);
32 crse(i,j,k) = a*b/(a+b);
33}
34
35AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
36void mlndlap_avgdown_coeff_y (int i, int j, int k, Array4<Real> const& crse,
37 Array4<Real const> const& fine) noexcept
38{
39 Real a = fine(2*i,2*j ,k) + fine(2*i+1,2*j ,k);
40 Real b = fine(2*i,2*j+1,k) + fine(2*i+1,2*j+1,k);
41 crse(i,j,k) = a*b/(a+b);
42}
43
44AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
45void mlndlap_semi_avgdown_coeff (int i, int j, int k, Array4<Real> const& crse,
46 Array4<Real const> const& fine, int idir) noexcept
47{
48 if (idir == 1) {
49 Real a = fine(2*i ,j,k);
50 Real b = fine(2*i+1,j,k);
51 crse(i,j,k) = Real(2.0)*a*b/(a+b);
52 } else {
53 Real a = fine(i,2*j ,k);
54 Real b = fine(i,2*j+1,k);
55 crse(i,j,k) = Real(2.0)*a*b/(a+b);
56 }
57}
58
59//
60// operator
61//
62
63AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64Real mlndlap_adotx_ha (int i, int j, int k, Array4<Real const> const& x,
65 Array4<Real const> const& sx, Array4<Real const> const& sy,
66 Array4<int const> const& msk, bool is_rz,
67 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
68{
69 if (msk(i,j,k)) {
70 return Real(0.0);
71 } else {
72 Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
73 Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
74 Real y = x(i-1,j-1,k)*(facx*sx(i-1,j-1,k)+facy*sy(i-1,j-1,k))
75 + x(i+1,j-1,k)*(facx*sx(i ,j-1,k)+facy*sy(i ,j-1,k))
76 + x(i-1,j+1,k)*(facx*sx(i-1,j ,k)+facy*sy(i-1,j ,k))
77 + x(i+1,j+1,k)*(facx*sx(i ,j ,k)+facy*sy(i ,j ,k))
78 + x(i-1,j,k)*(Real(2.0)*facx*(sx(i-1,j-1,k)+sx(i-1,j,k))
79 - facy*(sy(i-1,j-1,k)+sy(i-1,j,k)))
80 + x(i+1,j,k)*(Real(2.0)*facx*(sx(i ,j-1,k)+sx(i ,j,k))
81 - facy*(sy(i ,j-1,k)+sy(i ,j,k)))
82 + x(i,j-1,k)*( -facx*(sx(i-1,j-1,k)+sx(i,j-1,k))
83 +Real(2.0)*facy*(sy(i-1,j-1,k)+sy(i,j-1,k)))
84 + x(i,j+1,k)*( -facx*(sx(i-1,j ,k)+sx(i,j ,k))
85 +Real(2.0)*facy*(sy(i-1,j ,k)+sy(i,j ,k)))
86 + x(i,j,k)*(-Real(2.0))*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
87 +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
88 if (is_rz) {
89 Real fp = facy / static_cast<Real>(2*i+1);
90 Real fm = facy / static_cast<Real>(2*i-1);
91 y += (fm*sy(i-1,j ,k)-fp*sy(i,j ,k))*(x(i,j+1,k)-x(i,j,k))
92 + (fm*sy(i-1,j-1,k)-fp*sy(i,j-1,k))*(x(i,j-1,k)-x(i,j,k));
93 }
94 return y;
95 }
96}
97
98AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
99Real mlndlap_adotx_aa (int i, int j, int k, Array4<Real const> const& x,
100 Array4<Real const> const& sig, Array4<int const> const& msk,
101 bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
102{
103 if (msk(i,j,k)) {
104 return Real(0.0);
105 } else {
106 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
107 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
108 Real fxy = facx + facy;
109 Real f2xmy = Real(2.0)*facx - facy;
110 Real fmx2y = Real(2.0)*facy - facx;
111 Real y = x(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
112 + x(i+1,j-1,k)*fxy*sig(i ,j-1,k)
113 + x(i-1,j+1,k)*fxy*sig(i-1,j ,k)
114 + x(i+1,j+1,k)*fxy*sig(i ,j ,k)
115 + x(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
116 + x(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k))
117 + x(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
118 + x(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k))
119 + x(i,j,k)*(-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)
120 +sig(i-1,j,k)+sig(i,j,k));
121 if (is_rz) {
122 Real fp = facy / static_cast<Real>(2*i+1);
123 Real fm = facy / static_cast<Real>(2*i-1);
124 y += (fm*sig(i-1,j ,k)-fp*sig(i,j ,k))*(x(i,j+1,k)-x(i,j,k))
125 + (fm*sig(i-1,j-1,k)-fp*sig(i,j-1,k))*(x(i,j-1,k)-x(i,j,k));
126 }
127 return y;
128 }
129}
130
131AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
132Real mlndlap_adotx_c (int i, int j, int k, Array4<Real const> const& x,
133 Real sigma, Array4<int const> const& msk,
134 bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
135{
136 if (msk(i,j,k)) {
137 return Real(0.0);
138 } else {
139 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
140 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
141 Real fxy = facx + facy;
142 Real f2xmy = Real(2.0)*facx - facy;
143 Real fmx2y = Real(2.0)*facy - facx;
144 Real y = (x(i-1,j-1,k)*fxy
145 + x(i+1,j-1,k)*fxy
146 + x(i-1,j+1,k)*fxy
147 + x(i+1,j+1,k)*fxy
148 + x(i-1,j,k)*f2xmy*Real(2.)
149 + x(i+1,j,k)*f2xmy*Real(2.)
150 + x(i,j-1,k)*fmx2y*Real(2.)
151 + x(i,j+1,k)*fmx2y*Real(2.)
152 + x(i,j,k)*(-Real(2.0))*fxy*Real(4.));
153 if (is_rz) {
154 Real fp = facy / static_cast<Real>(2*i+1);
155 Real fm = facy / static_cast<Real>(2*i-1);
156 y += ((fm-fp)*(x(i,j+1,k)-x(i,j,k))
157 + (fm-fp)*(x(i,j-1,k)-x(i,j,k)));
158 }
159 return y * sigma;
160 }
161}
162
163AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
164void mlndlap_normalize_ha (int i, int j, int k, Array4<Real> const& x, Array4<Real const> const& sx,
165 Array4<Real const> const& sy, Array4<int const> const& msk,
166 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
167{
168 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
169 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
170
171 if (!msk(i,j,k)) {
172 x(i,j,k) /= (-Real(2.0))*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
173 +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
174 }
175}
176
177AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
178void mlndlap_normalize_aa (int i, int j, int k, Array4<Real> const& x, Array4<Real const> const& sig,
179 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
180{
181 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
182 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
183 Real fxy = facx + facy;
184
185 if (!msk(i,j,k)) {
186 x(i,j,k) /= (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
187 }
188}
189
190AMREX_GPU_DEVICE AMREX_FORCE_INLINE
191void mlndlap_jacobi_ha (int i, int j, int k, Array4<Real> const& sol, Real Ax,
192 Array4<Real const> const& rhs, Array4<Real const> const& sx,
193 Array4<Real const> const& sy, Array4<int const> const& msk,
194 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
195{
196 Real facx = -Real(2.0/6.0)*dxinv[0]*dxinv[0];
197 Real facy = -Real(2.0/6.0)*dxinv[1]*dxinv[1];
198
199 if (msk(i,j,k)) {
200 sol(i,j,k) = Real(0.0);
201 } else {
202 sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
203 / (facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
204 + facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
205 }
206}
207
208inline
209void mlndlap_jacobi_ha (Box const& bx, Array4<Real> const& sol, Array4<Real const> const& Ax,
210 Array4<Real const> const& rhs, Array4<Real const> const& sx,
211 Array4<Real const> const& sy, Array4<int const> const& msk,
212 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
213{
214 Real facx = -Real(2.0/6.0)*dxinv[0]*dxinv[0];
215 Real facy = -Real(2.0/6.0)*dxinv[1]*dxinv[1];
216
217 amrex::LoopConcurrentOnCpu(bx, [&] (int i, int j, int k) noexcept
218 {
219 if (msk(i,j,k)) {
220 sol(i,j,k) = Real(0.0);
221 } else {
222 sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax(i,j,k))
223 / (facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
224 + facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
225 }
226 });
227}
228
229AMREX_GPU_DEVICE AMREX_FORCE_INLINE
230void mlndlap_jacobi_aa (int i, int j, int k, Array4<Real> const& sol, Real Ax,
231 Array4<Real const> const& rhs, Array4<Real const> const& sig,
232 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
233{
234 Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
235
236 if (msk(i,j,k)) {
237 sol(i,j,k) = Real(0.0);
238 } else {
239 sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
240 / (fac*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k)));
241 }
242}
243
244AMREX_GPU_DEVICE AMREX_FORCE_INLINE
245void mlndlap_jacobi_c (int i, int j, int k, Array4<Real> const& sol, Real Ax,
246 Array4<Real const> const& rhs, Real sig,
247 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
248{
249 Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
250
251 if (msk(i,j,k)) {
252 sol(i,j,k) = Real(0.0);
253 } else {
254 sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax)
255 / (fac*Real(4.)*sig);
256 }
257}
258
259inline
260void mlndlap_jacobi_aa (Box const& bx, Array4<Real> const& sol, Array4<Real const> const& Ax,
261 Array4<Real const> const& rhs, Array4<Real const> const& sig,
262 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
263{
264 Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
265
266 amrex::LoopConcurrentOnCpu(bx, [&] (int i, int j, int k) noexcept
267 {
268 if (msk(i,j,k)) {
269 sol(i,j,k) = Real(0.0);
270 } else {
271 sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax(i,j,k))
272 / (fac*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k)));
273 }
274 });
275}
276
277inline
278void mlndlap_jacobi_c (Box const& bx, Array4<Real> const& sol, Array4<Real const> const& Ax,
279 Array4<Real const> const& rhs, Real sig,
280 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
281{
282 Real fac = -Real(2.0/6.0)*(dxinv[0]*dxinv[0] + dxinv[1]*dxinv[1]);
283
284 amrex::LoopConcurrentOnCpu(bx, [&] (int i, int j, int k) noexcept
285 {
286 if (msk(i,j,k)) {
287 sol(i,j,k) = Real(0.0);
288 } else {
289 sol(i,j,k) += Real(2.0/3.0) * (rhs(i,j,k) - Ax(i,j,k))
290 / (fac*Real(4.)*sig);
291 }
292 });
293}
294
295inline
296void mlndlap_gauss_seidel_ha (Box const& bx, Array4<Real> const& sol,
297 Array4<Real const> const& rhs, Array4<Real const> const& sx,
298 Array4<Real const> const& sy, Array4<int const> const& msk,
299 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
300 bool is_rz) noexcept
301{
302 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
303 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
304
305 amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
306 {
307 if (msk(i,j,k)) {
308 sol(i,j,k) = Real(0.0);
309 } else {
310 Real s0 = Real(-2.0)*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
311 +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
312
313 Real Ax = sol(i-1,j-1,k)*(facx*sx(i-1,j-1,k)+facy*sy(i-1,j-1,k))
314 + sol(i+1,j-1,k)*(facx*sx(i ,j-1,k)+facy*sy(i ,j-1,k))
315 + sol(i-1,j+1,k)*(facx*sx(i-1,j ,k)+facy*sy(i-1,j ,k))
316 + sol(i+1,j+1,k)*(facx*sx(i ,j ,k)+facy*sy(i ,j ,k))
317 + sol(i-1,j,k)*(Real(2.0)*facx*(sx(i-1,j-1,k)+sx(i-1,j,k))
318 - facy*(sy(i-1,j-1,k)+sy(i-1,j,k)))
319 + sol(i+1,j,k)*(Real(2.0)*facx*(sx(i ,j-1,k)+sx(i ,j,k))
320 - facy*(sy(i ,j-1,k)+sy(i ,j,k)))
321 + sol(i,j-1,k)*( -facx*(sx(i-1,j-1,k)+sx(i,j-1,k))
322 +Real(2.0)*facy*(sy(i-1,j-1,k)+sy(i,j-1,k)))
323 + sol(i,j+1,k)*( -facx*(sx(i-1,j ,k)+sx(i,j ,k))
324 +Real(2.0)*facy*(sy(i-1,j ,k)+sy(i,j ,k)))
325 + sol(i,j,k)*s0;
326
327 if (is_rz) {
328 Real fp = facy / static_cast<Real>(2*i+1);
329 Real fm = facy / static_cast<Real>(2*i-1);
330 Real frzlo = fm*sy(i-1,j-1,k)-fp*sy(i,j-1,k);
331 Real frzhi = fm*sy(i-1,j ,k)-fp*sy(i,j ,k);
332 s0 += - frzhi - frzlo;
333 Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
334 + frzlo*(sol(i,j-1,k)-sol(i,j,k));
335 }
336
337 sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
338 }
339 });
340}
341
342inline
343void mlndlap_gauss_seidel_aa (Box const& bx, Array4<Real> const& sol,
344 Array4<Real const> const& rhs, Array4<Real const> const& sig,
345 Array4<int const> const& msk,
346 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
347 bool is_rz) noexcept
348{
349 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
350 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
351 Real fxy = facx + facy;
352 Real f2xmy = Real(2.0)*facx - facy;
353 Real fmx2y = Real(2.0)*facy - facx;
354
355 amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
356 {
357 if (msk(i,j,k)) {
358 sol(i,j,k) = Real(0.0);
359 } else {
360 Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
361 Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
362 + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
363 + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
364 + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
365 + sol(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
366 + sol(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k))
367 + sol(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
368 + sol(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k))
369 + sol(i,j,k)*s0;
370
371 if (is_rz) {
372 Real fp = facy / static_cast<Real>(2*i+1);
373 Real fm = facy / static_cast<Real>(2*i-1);
374 Real frzlo = fm*sig(i-1,j-1,k)-fp*sig(i,j-1,k);
375 Real frzhi = fm*sig(i-1,j ,k)-fp*sig(i,j ,k);
376 s0 += - frzhi - frzlo;
377 Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
378 + frzlo*(sol(i,j-1,k)-sol(i,j,k));
379 }
380
381 sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
382 }
383 });
384}
385
386inline
387void mlndlap_gauss_seidel_c (Box const& bx, Array4<Real> const& sol,
388 Array4<Real const> const& rhs, Real sig,
389 Array4<int const> const& msk,
390 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
391 bool is_rz) noexcept
392{
393 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
394 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
395 Real fxy = facx + facy;
396 Real f2xmy = Real(2.0)*facx - facy;
397 Real fmx2y = Real(2.0)*facy - facx;
398
399 amrex::LoopOnCpu(bx, [&] (int i, int j, int k) noexcept
400 {
401 if (msk(i,j,k)) {
402 sol(i,j,k) = Real(0.0);
403 } else {
404 Real s0 = (-Real(2.0))*fxy*Real(4.);
405 Real Ax = sol(i-1,j-1,k)*fxy
406 + sol(i+1,j-1,k)*fxy
407 + sol(i-1,j+1,k)*fxy
408 + sol(i+1,j+1,k)*fxy
409 + sol(i-1,j,k)*f2xmy*Real(2.)
410 + sol(i+1,j,k)*f2xmy*Real(2.)
411 + sol(i,j-1,k)*fmx2y*Real(2.)
412 + sol(i,j+1,k)*fmx2y*Real(2.)
413 + sol(i,j,k)*s0;
414
415 if (is_rz) {
416 Real fp = facy / static_cast<Real>(2*i+1);
417 Real fm = facy / static_cast<Real>(2*i-1);
418 Real frzlo = fm-fp;
419 Real frzhi = fm-fp;
420 s0 += - frzhi - frzlo;
421 Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
422 + frzlo*(sol(i,j-1,k)-sol(i,j,k));
423 }
424
425 sol(i,j,k) += (rhs(i,j,k) - Ax*sig) / (s0*sig);
426 }
427 });
428}
429
430AMREX_FORCE_INLINE
431void tridiagonal_solve (Array1D<Real,0,31>& a_ls, Array1D<Real,0,31>& b_ls, Array1D<Real,0,31>& c_ls,
432 Array1D<Real,0,31>& r_ls, Array1D<Real,0,31>& u_ls, Array1D<Real,0,31>& gam,
433 int ilen ) noexcept
434{
435 Real bet = b_ls(0);
436 u_ls(0) = r_ls(0) / bet;
437
438 for (int i = 1; i <= ilen - 1; i++) {
439 gam(i) = c_ls(i-1) / bet;
440 bet = b_ls(i) - a_ls(i)*gam(i);
441 if (bet == 0) { amrex::Abort(">>>TRIDIAG FAILED"); }
442 u_ls(i) = (r_ls(i)-a_ls(i)*u_ls(i-1)) / bet;
443 }
444 for (int i = ilen-2; i >= 0; i--) {
445 u_ls(i) = u_ls(i) - gam(i+1)*u_ls(i+1);
446 }
447}
448
449inline
450void mlndlap_gauss_seidel_with_line_solve_aa (Box const& bx, Array4<Real> const& sol,
451 Array4<Real const> const& rhs, Array4<Real const> const& sig,
452 Array4<int const> const& msk,
453 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
454 bool is_rz) noexcept
455{
456 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
457 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
458 Real fxy = facx + facy;
459 Real f2xmy = Real(2.0)*facx - facy;
460 Real fmx2y = Real(2.0)*facy - facx;
461
462 if (is_rz) {
463 amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa is not implemented in r-z 2D ");
464 }
465
466 const auto lo = amrex::lbound(bx);
467 const auto hi = amrex::ubound(bx);
468
469 int idir = -1;
470 int ilen = 33;
471 int k = 0;
472 if (dxinv[0] <= dxinv[1]) {
473 idir = 1;
474 ilen = hi.y - lo.y + 1;
475 }
476 if (dxinv[1] <= dxinv[0]) {
477 idir = 0;
478 ilen = hi.x - lo.x + 1;
479 }
480
481 if (ilen > 32) {
482 amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa is hard-wired to be no longer than 32");
483 }
484
485 Array1D<Real,0,31> a_ls,b_ls,c_ls,u_ls,r_ls,gam;
486
487 if (idir == 1) {
488 for (int i = lo.x; i <= hi.x; ++i)
489 {
490 for (int j = lo.y; j <= hi.y; ++j)
491 {
492 if (msk(i,j,k)) {
493 a_ls(j-lo.y) = 0.;
494 b_ls(j-lo.y) = 1.;
495 c_ls(j-lo.y) = 0.;
496 u_ls(j-lo.y) = 0.;
497 r_ls(j-lo.y) = 0.;
498 }
499 else
500 {
501 Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
502
503 Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
504 + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
505 + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
506 + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
507 + sol(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
508 + sol(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k));
509
510 a_ls(j-lo.y) = fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k));
511 b_ls(j-lo.y) = s0;
512 c_ls(j-lo.y) = fmx2y*(sig(i-1,j ,k)+sig(i,j ,k));
513 u_ls(j-lo.y) = 0.;
514 r_ls(j-lo.y) = rhs(i,j,k) - Ax;
515 }
516 }
517 tridiagonal_solve(a_ls, b_ls, c_ls, r_ls, u_ls, gam, ilen);
518
519 for (int j = lo.y; j <= hi.y; ++j)
520 {
521 sol(i,j,k) = u_ls(j-lo.y);
522 }
523 }
524 } else if (idir == 0) {
525 for (int j = lo.y ;j <= hi.y; ++j)
526 {
527 for (int i = lo.x; i <= hi.x; ++i)
528 {
529 if (msk(i,j,k)) {
530 a_ls(i-lo.x) = 0.;
531 b_ls(i-lo.x) = 1.;
532 c_ls(i-lo.x) = 0.;
533 u_ls(i-lo.x) = 0.;
534 r_ls(i-lo.x) = 0.;
535 }
536 else
537 {
538 Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
539
540 Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
541 + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
542 + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
543 + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
544 + sol(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
545 + sol(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k));
546
547 a_ls(i-lo.x) = f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k));
548 b_ls(i-lo.x) = s0;
549 c_ls(i-lo.x) = f2xmy*(sig(i ,j-1,k)+sig(i ,j,k));
550 u_ls(i-lo.x) = 0.;
551 r_ls(i-lo.x) = rhs(i,j,k) - Ax;
552
553 }
554 }
555 tridiagonal_solve(a_ls, b_ls, c_ls, r_ls, u_ls, gam, ilen);
556
557 for (int i = lo.x; i <= hi.x; ++i)
558 {
559 sol(i,j,k) = u_ls(i-lo.x);
560 }
561 }
562 } else {
563 amrex::Abort("mlndlap_gauss_seidel_with_line_solve_aa is wrong direction.");
564 }
565
566}
567
568//
569// interpolation
570//
571
572 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
573 Real aa_interp_line_x (Array4<Real const> const& crse, Array4<Real const> const& sig,
574 int i, int j, int ic, int jc) noexcept
575 {
576 Real w1 = sig(i-1,j-1,0) + sig(i-1,j,0);
577 Real w2 = sig(i ,j-1,0) + sig(i ,j,0);
578 return (w1*crse(ic,jc,0)+w2*crse(ic+1,jc,0))/(w1+w2);
579 }
580
581 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
582 Real aa_interp_line_y (Array4<Real const> const& crse, Array4<Real const> const& sig,
583 int i, int j, int ic, int jc) noexcept
584 {
585 Real w1 = sig(i-1,j-1,0) + sig(i,j-1,0);
586 Real w2 = sig(i-1,j ,0) + sig(i,j ,0);
587 return (w1*crse(ic,jc,0)+w2*crse(ic,jc+1,0))/(w1+w2);
588 }
589
590 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
591 Real aa_interp_face_xy (Array4<Real const> const& crse, Array4<Real const> const& sig,
592 int i, int j, int ic, int jc) noexcept
593 {
594 Real w1 = sig(i-1,j-1,0) + sig(i-1,j,0);
595 Real w2 = sig(i ,j-1,0) + sig(i ,j,0);
596 Real w3 = sig(i-1,j-1,0) + sig(i,j-1,0);
597 Real w4 = sig(i-1,j ,0) + sig(i,j ,0);
598 return (w1 * aa_interp_line_y(crse,sig,i-1,j ,ic ,jc ) +
599 w2 * aa_interp_line_y(crse,sig,i+1,j ,ic+1,jc ) +
600 w3 * aa_interp_line_x(crse,sig,i ,j-1,ic ,jc ) +
601 w4 * aa_interp_line_x(crse,sig,i ,j+1,ic ,jc+1)) / (w1+w2+w3+w4);
602 }
603
604AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
605void mlndlap_interpadd_c (int i, int j, int, Array4<Real> const& fine,
606 Array4<Real const> const& crse,
607 Array4<int const> const& msk) noexcept
608{
609 if (!msk(i,j,0)) {
610 int ic = amrex::coarsen(i,2);
611 int jc = amrex::coarsen(j,2);
612 bool i_is_odd = (ic*2 != i);
613 bool j_is_odd = (jc*2 != j);
614 if (i_is_odd && j_is_odd) {
615 // Node on a X-Y face
616 fine(i,j,0) += Real(0.25) * (crse(ic ,jc ,0) +
617 crse(ic+1,jc ,0) +
618 crse(ic ,jc+1,0) +
619 crse(ic+1,jc+1,0));
620 } else if (i_is_odd) {
621 // Node on X line
622 fine(i,j,0) += Real(0.5) * (crse(ic,jc,0)+crse(ic+1,jc,0));
623 } else if (j_is_odd) {
624 // Node on Y line
625 fine(i,j,0) += Real(0.5) * (crse(ic,jc,0)+crse(ic,jc+1,0));
626 } else {
627 // Node coincident with coarse node
628 fine(i,j,0) += crse(ic,jc,0);
629 }
630 }
631}
632
633AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
634void mlndlap_interpadd_aa (int i, int j, int, Array4<Real> const& fine,
635 Array4<Real const> const& crse, Array4<Real const> const& sig,
636 Array4<int const> const& msk) noexcept
637{
638 if (!msk(i,j,0)) {
639 int ic = amrex::coarsen(i,2);
640 int jc = amrex::coarsen(j,2);
641 bool i_is_odd = (ic*2 != i);
642 bool j_is_odd = (jc*2 != j);
643 if (i_is_odd && j_is_odd) {
644 // Node on a X-Y face
645 fine(i,j,0) += aa_interp_face_xy(crse,sig,i,j,ic,jc);
646 } else if (i_is_odd) {
647 // Node on X line
648 fine(i,j,0) += aa_interp_line_x(crse,sig,i,j,ic,jc);
649 } else if (j_is_odd) {
650 // Node on Y line
651 fine(i,j,0) += aa_interp_line_y(crse,sig,i,j,ic,jc);
652 } else {
653 // Node coincident with coarse node
654 fine(i,j,0) += crse(ic,jc,0);
655 }
656 }
657}
658
659AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
660void mlndlap_semi_interpadd_aa (int i, int j, int, Array4<Real> const& fine,
661 Array4<Real const> const& crse, Array4<Real const> const& sig,
662 Array4<int const> const& msk, int idir) noexcept
663{
664 if (idir == 1) {
665 if (!msk(i,j,0)) {
666 int ic = amrex::coarsen(i,2);
667 int jc = j;
668 bool i_is_odd = (ic*2 != i);
669 if (i_is_odd) {
670 // Node on X line
671 fine(i,j,0) += aa_interp_line_x(crse,sig,i,j,ic,jc);
672 } else {
673 //Node coincident with coarse node
674 fine(i,j,0) += crse(ic,jc,0);
675 }
676 }
677 } else if (idir == 0 ) {
678 if (!msk(i,j,0)) {
679 int ic = i;
680 int jc = amrex::coarsen(j,2);
681 bool j_is_odd = (ic*2 != i);
682 if (j_is_odd) {
683 // Node on Y line
684 fine(i,j,0) += aa_interp_line_y(crse,sig,i,j,ic,jc);
685 } else {
686 //Node coincident with coarse node
687 fine(i,j,0) += crse(ic,jc,0);
688 }
689 }
690 }
691}
692
693 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
694 Real ha_interp_face_xy (Array4<Real const> const& crse,
695 Array4<Real const> const& sigx, Array4<Real const> const& sigy,
696 int i, int j, int ic, int jc) noexcept
697 {
698 Real w1 = sigx(i-1,j-1,0) + sigx(i-1,j,0);
699 Real w2 = sigx(i ,j-1,0) + sigx(i ,j,0);
700 Real w3 = sigy(i-1,j-1,0) + sigy(i,j-1,0);
701 Real w4 = sigy(i-1,j ,0) + sigy(i,j ,0);
702 return (w1 * aa_interp_line_y(crse,sigy,i-1,j ,ic ,jc ) +
703 w2 * aa_interp_line_y(crse,sigy,i+1,j ,ic+1,jc ) +
704 w3 * aa_interp_line_x(crse,sigx,i ,j-1,ic ,jc ) +
705 w4 * aa_interp_line_x(crse,sigx,i ,j+1,ic ,jc+1)) / (w1+w2+w3+w4);
706 }
707
708AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
709void mlndlap_interpadd_ha (int i, int j, int,
710 Array4<Real> const& fine, Array4<Real const> const& crse,
711 Array4<Real const> const& sigx, Array4<Real const> const& sigy,
712 Array4<int const> const& msk) noexcept
713{
714 if (!msk(i,j,0)) {
715 int ic = amrex::coarsen(i,2);
716 int jc = amrex::coarsen(j,2);
717 bool i_is_odd = (ic*2 != i);
718 bool j_is_odd = (jc*2 != j);
719 if (i_is_odd && j_is_odd) {
720 // Node on a X-Y face
721 fine(i,j,0) += ha_interp_face_xy(crse,sigx,sigy,i,j,ic,jc);
722 } else if (i_is_odd) {
723 // Node on X line
724 fine(i,j,0) += aa_interp_line_x(crse,sigx,i,j,ic,jc);
725 } else if (j_is_odd) {
726 // Node on Y line
727 fine(i,j,0) += aa_interp_line_y(crse,sigy,i,j,ic,jc);
728 } else {
729 // Node coincident with coarse node
730 fine(i,j,0) += crse(ic,jc,0);
731 }
732 }
733}
734
735//
736// rhs & u
737//
738
739AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
740void mlndlap_divu (int i, int j, int k, Array4<Real> const& rhs, Array4<Real const> const& vel,
741 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
742 Box const& nodal_domain,
743 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
744 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
745 bool is_rz) noexcept
746{
747 Real facx = Real(0.5)*dxinv[0];
748 Real facy = Real(0.5)*dxinv[1];
749
750 const auto domlo = amrex::lbound(nodal_domain);
751 const auto domhi = amrex::ubound(nodal_domain);
752
753 if (msk(i,j,k)) {
754 rhs(i,j,k) = Real(0.0);
755 } else {
756
757 Real zero_ilo = Real(1.0);
758 Real zero_ihi = Real(1.0);
759 Real zero_jlo = Real(1.0);
760 Real zero_jhi = Real(1.0);
761
762 // The nodal divergence operator should not see the tangential velocity
763 // at an inflow face
764 if ((bclo[0] == LinOpBCType::Neumann || bclo[0] == LinOpBCType::inflow)
765 && i == domlo.x)
766 {
767 zero_ilo = Real(0.0);
768 }
769 if ((bchi[0] == LinOpBCType::Neumann || bchi[0] == LinOpBCType::inflow)
770 && i == domhi.x)
771 {
772 zero_ihi = Real(0.0);
773 }
774 if ((bclo[1] == LinOpBCType::Neumann || bclo[1] == LinOpBCType::inflow)
775 && j == domlo.y)
776 {
777 zero_jlo = Real(0.0);
778 }
779 if ((bchi[1] == LinOpBCType::Neumann || bchi[1] == LinOpBCType::inflow)
780 && j == domhi.y)
781 {
782 zero_jhi = Real(0.0);
783 }
784
785 rhs(i,j,k) = facx*(-vel(i-1,j-1,k,0)*zero_jlo + vel(i,j-1,k,0)*zero_jlo
786 -vel(i-1,j ,k,0)*zero_jhi + vel(i,j ,k,0)*zero_jhi)
787 + facy*(-vel(i-1,j-1,k,1)*zero_ilo - vel(i,j-1,k,1)*zero_ihi
788 +vel(i-1,j ,k,1)*zero_ilo + vel(i,j ,k,1)*zero_ihi);
789 if (is_rz) {
790 // Here we assume we can't have inflow in the radial direction
791 Real fm = facy / static_cast<Real>(6*i-3);
792 Real fp = facy / static_cast<Real>(6*i+3);
793 rhs(i,j,k) += fm*(vel(i-1,j,k,1)-vel(i-1,j-1,k,1))
794 - fp*(vel(i ,j,k,1)-vel(i ,j-1,k,1));
795 }
796 }
797}
798
799AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
800Real mlndlap_rhcc (int i, int j, int k, Array4<Real const> const& rhcc,
801 Array4<int const> const& msk) noexcept
802{
803 Real r;
804 if (msk(i,j,k)) {
805 r = Real(0.0);
806 } else {
807 r = Real(0.25) * (rhcc(i-1,j-1,k)+rhcc(i,j-1,k)+rhcc(i-1,j,k)+rhcc(i,j,k));
808 }
809 return r;
810}
811
812AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
813void mlndlap_mknewu (int i, int j, int k, Array4<Real> const& u, Array4<Real const> const& p,
814 Array4<Real const> const& sig, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
815 bool is_rz) noexcept
816{
817 Real facx = Real(0.5)*dxinv[0];
818 Real facy = Real(0.5)*dxinv[1];
819 u(i,j,k,0) -= sig(i,j,k)*facx*(-p(i,j,k)+p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
820 u(i,j,k,1) -= sig(i,j,k)*facy*(-p(i,j,k)-p(i+1,j,k)+p(i,j+1,k)+p(i+1,j+1,k));
821 if (is_rz) {
822 Real frz = sig(i,j,k)*facy / static_cast<Real>(6*i+3);
823 u(i,j,k,1) += frz*(p(i,j,k)-p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
824 }
825}
826
827AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
828void mlndlap_mknewu_c (int i, int j, int k, Array4<Real> const& u, Array4<Real const> const& p,
829 Real sig, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
830 bool is_rz) noexcept
831{
832 Real facx = Real(0.5)*dxinv[0];
833 Real facy = Real(0.5)*dxinv[1];
834 u(i,j,k,0) -= sig*facx*(-p(i,j,k)+p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
835 u(i,j,k,1) -= sig*facy*(-p(i,j,k)-p(i+1,j,k)+p(i,j+1,k)+p(i+1,j+1,k));
836 if (is_rz) {
837 Real frz = sig*facy / static_cast<Real>(6*i+3);
838 u(i,j,k,1) += frz*(p(i,j,k)-p(i+1,j,k)-p(i,j+1,k)+p(i+1,j+1,k));
839 }
840}
841
842 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
843 Real mlndlap_sum_Df (int ii, int jj, Real facx, Real facy,
844 Array4<Real const> const& vel, Box const& velbx, bool is_rz) noexcept
845 {
846 Real fm = is_rz ? facy / static_cast<Real>(6*ii-3) : Real(0.0);
847 Real fp = is_rz ? facy / static_cast<Real>(6*ii+3) : Real(0.0);
848
849 Real Df = Real(0.0);
850 if (velbx.contains(ii-1,jj-1,0)) {
851 Df += -facx*vel(ii-1,jj-1,0,0) - (facy+fm)*vel(ii-1,jj-1,0,1);
852 }
853 if (velbx.contains(ii,jj-1,0)) {
854 Df += facx*vel(ii,jj-1,0,0) - (facy-fp)*vel(ii,jj-1,0,1);
855 }
856 if (velbx.contains(ii-1,jj,0)) {
857 Df += -facx*vel(ii-1,jj,0,0) + (facy+fm)*vel(ii-1,jj,0,1);
858 }
859 if (velbx.contains(ii,jj,0)) {
860 Df += facx*vel(ii,jj,0,0) + (facy-fp)*vel(ii,jj,0,1);
861 }
862 return Df;
863 }
864
865template <int rr>
866AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
867void mlndlap_divu_fine_contrib (int i, int j, int /*k*/, Box const& fvbx, Box const& velbx,
868 Array4<Real> const& rhs, Array4<Real const> const& vel,
869 Array4<Real const> const& frhs, Array4<int const> const& msk,
870 GpuArray<Real,AMREX_SPACEDIM> const& dxinv, bool is_rz) noexcept
871{
872 const int ii = rr*i;
873 const int jj = rr*j;
874 if (msk(ii,jj,0)) {
875 const Real facx = Real(0.5)*dxinv[0];
876 const Real facy = Real(0.5)*dxinv[1];
877
878 Real Df = Real(0.0);
879
880 const int ilo = amrex::max(ii-rr+1, fvbx.smallEnd(0));
881 const int ihi = amrex::min(ii+rr-1, fvbx.bigEnd (0));
882 const int jlo = amrex::max(jj-rr+1, fvbx.smallEnd(1));
883 const int jhi = amrex::min(jj+rr-1, fvbx.bigEnd (1));
884
885 for (int joff = jlo; joff <= jhi; ++joff) {
886 for (int ioff = ilo; ioff <= ihi; ++ioff) {
887 Real scale = static_cast<Real>((rr-std::abs(ii-ioff)) *
888 (rr-std::abs(jj-joff)));
889 if (fvbx.strictly_contains(ioff,joff,0)) {
890 Df += scale * frhs(ioff,joff,0);
891 } else {
892 Df += scale * mlndlap_sum_Df(ioff, joff, facx, facy, vel, velbx, is_rz);
893 }
894 }}
895
896 rhs(i,j,0) = Df * (Real(1.0)/static_cast<Real>(rr*rr*rr*rr));
897 } else {
898 rhs(i,j,0) = Real(0.0);
899 }
900}
901
902template <int rr>
903AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
904void mlndlap_rhcc_fine_contrib (int i, int j, int, Box const& ccbx,
905 Array4<Real> const& rhs, Array4<Real const> const& cc,
906 Array4<int const> const& msk) noexcept
907{
908 const int ii = rr*i;
909 const int jj = rr*j;
910 if (msk(ii,jj,0)) {
911 Real tmp = Real(0.0);
912
913 const int ilo = amrex::max(ii-rr , ccbx.smallEnd(0));
914 const int ihi = amrex::min(ii+rr-1, ccbx.bigEnd (0));
915 const int jlo = amrex::max(jj-rr , ccbx.smallEnd(1));
916 const int jhi = amrex::min(jj+rr-1, ccbx.bigEnd (1));
917
918 for (int joff = jlo; joff <= jhi; ++joff) {
919 for (int ioff = ilo; ioff <= ihi; ++ioff) {
920 Real scale = (static_cast<Real>(rr)-std::abs(static_cast<Real>(ioff-ii)+Real(0.5)))
921 * (static_cast<Real>(rr)-std::abs(static_cast<Real>(joff-jj)+Real(0.5)));
922 tmp += cc(ioff,joff,0) * scale;
923 }}
924
925 rhs(i,j,0) += tmp * (Real(1.0)/Real(rr*rr*rr*rr));
926 }
927}
928
929 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
930 Real neumann_scale (int i, int j, Box const& nddom,
931 GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bclo,
932 GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bchi) noexcept
933 {
934 Real val = Real(1.0);
935
936 const auto ndlo = amrex::lbound(nddom);
937 const auto ndhi = amrex::ubound(nddom);
938
939 if ((i == ndlo.x && ( bclo[0] == LinOpBCType::Neumann ||
940 bclo[0] == LinOpBCType::inflow)) ||
941 (i == ndhi.x && ( bchi[0] == LinOpBCType::Neumann ||
942 bchi[0] == LinOpBCType::inflow))) {
943 val *= Real(2.0);
944 }
945
946 if ((j == ndlo.y && ( bclo[1] == LinOpBCType::Neumann ||
947 bclo[1] == LinOpBCType::inflow)) ||
948 (j == ndhi.y && ( bchi[1] == LinOpBCType::Neumann ||
949 bchi[1] == LinOpBCType::inflow))) {
950 val *= Real(2.0);
951 }
952
953 return val;
954 }
955
956AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
957void mlndlap_divu_cf_contrib (int i, int j, int, Array4<Real> const& rhs,
958 Array4<Real const> const& vel, Array4<Real const> const& fc,
959 Array4<Real const> const& rhcc, Array4<int const> const& dmsk,
960 Array4<int const> const& ndmsk, Array4<int const> const& ccmsk,
961 bool is_rz, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
962 Box const& ccdom_p, Box const& veldom, Box const& nddom,
963 GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bclo,
964 GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bchi) noexcept
965{
966 using namespace nodelap_detail;
967
968 if (!dmsk(i,j,0) && ndmsk(i,j,0) == crse_fine_node) {
969 Real facx = Real(0.5) * dxinv[0];
970 Real facy = Real(0.5) * dxinv[1];
971 Real tmp = fc(i,j,0);
972
973 Real fm = is_rz ? facy / static_cast<Real>(6*i-3) : Real(0.0);
974 Real fp = is_rz ? facy / static_cast<Real>(6*i+3) : Real(0.0);
975
976 // Where there is inflow, veldom there is bigger than ccdom_p by one cell.
977 // ccdom_p is cc domain grown at periodic boundaries.
978
979 if (ccmsk(i-1,j-1,0) == crse_cell && veldom.contains(i-1,j-1,0)) {
980 tmp += -facx*vel(i-1,j-1,0,0) - (facy+fm)*vel(i-1,j-1,0,1);
981 if (rhcc && ccdom_p.contains(i-1,j-1,0)) {
982 tmp += Real(0.25) * rhcc(i-1,j-1,0);
983 }
984 }
985
986 if (ccmsk(i,j-1,0) == crse_cell && veldom.contains(i,j-1,0)) {
987 tmp += facx*vel(i,j-1,0,0) - (facy-fp)*vel(i,j-1,0,1);
988 if (rhcc && ccdom_p.contains(i,j-1,0)) {
989 tmp += Real(0.25) * rhcc(i,j-1,0);
990 }
991 }
992
993 if (ccmsk(i-1,j,0) == crse_cell && veldom.contains(i-1,j,0)) {
994 tmp += -facx*vel(i-1,j,0,0) + (facy+fm)*vel(i-1,j,0,1);
995 if (rhcc && ccdom_p.contains(i-1,j,0)) {
996 tmp += Real(0.25) * rhcc(i-1,j,0);
997 }
998 }
999
1000 if (ccmsk(i,j,0) == crse_cell && veldom.contains(i,j,0)) {
1001 tmp += facx*vel(i,j,0,0) + (facy-fp)*vel(i,j,0,1);
1002 if (rhcc && ccdom_p.contains(i,j,0)) {
1003 tmp += Real(0.25) * rhcc(i,j,0);
1004 }
1005 }
1006
1007 rhs(i,j,0) = tmp * neumann_scale(i, j, nddom, bclo, bchi);
1008 }
1009}
1010
1011//
1012// residual
1013//
1014
1015AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1016void mlndlap_crse_resid (int i, int j, int k, Array4<Real> const& resid,
1017 Array4<Real const> const& rhs, Array4<int const> const& msk,
1018 Box const& nddom, GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bclo,
1019 GpuArray<LinOpBCType,AMREX_SPACEDIM> const& bchi,
1020 bool neumann_doubling) noexcept
1021{
1022 if ( msk(i-1,j-1,k ) == 0 ||
1023 msk(i ,j-1,k ) == 0 ||
1024 msk(i-1,j ,k ) == 0 ||
1025 msk(i ,j ,k ) == 0 )
1026 {
1027 Real fac = Real(1.0);
1028 if (neumann_doubling) {
1029 const auto ndlo = amrex::lbound(nddom);
1030 const auto ndhi = amrex::ubound(nddom);
1031 if ((i == ndlo.x && ( bclo[0] == LinOpBCType::Neumann ||
1032 bclo[0] == LinOpBCType::inflow)) ||
1033 (i == ndhi.x && ( bchi[0] == LinOpBCType::Neumann ||
1034 bchi[0] == LinOpBCType::inflow))) {
1035 fac *= Real(2.);
1036 }
1037 if ((j == ndlo.y && ( bclo[1] == LinOpBCType::Neumann ||
1038 bclo[1] == LinOpBCType::inflow)) ||
1039 (j == ndhi.y && ( bchi[1] == LinOpBCType::Neumann ||
1040 bchi[1] == LinOpBCType::inflow))) {
1041 fac *= Real(2.);
1042 }
1043 }
1044
1045 resid(i,j,k) = (rhs(i,j,k) - resid(i,j,k)) * fac;
1046 } else {
1047 resid(i,j,k) = Real(0.);
1048 }
1049}
1050
1051//
1052// sync residual
1053//
1054
1055 template <typename P, typename S>
1056 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1057 Real mlndlap_sum_Ax (P const& pred, S const& sig, int i, int j, Real facx, Real facy,
1058 Array4<Real const> const& phi, bool is_rz) noexcept
1059 {
1060 Real Ax = Real(0.0);
1061 Real fp = Real(0.0), fm = Real(0.0);
1062 if (is_rz) {
1063 fp = facy / static_cast<Real>(2*i+1);
1064 fm = facy / static_cast<Real>(2*i-1);
1065 }
1066 if (pred(i-1,j-1)) {
1067 Ax += sig(i-1,j-1,0)*(facx*(Real(2.)*(phi(i-1,j ,0)-phi(i ,j ,0))
1068 + (phi(i-1,j-1,0)-phi(i ,j-1,0)))
1069 + facy*(Real(2.)*(phi(i ,j-1,0)-phi(i ,j ,0))
1070 + (phi(i-1,j-1,0)-phi(i-1,j ,0)))
1071 + fm * (phi(i ,j-1,0)-phi(i ,j ,0)));
1072 }
1073 if (pred(i,j-1)) {
1074 Ax += sig(i,j-1,0)*(facx*(Real(2.)*(phi(i+1,j ,0)-phi(i ,j ,0))
1075 + (phi(i+1,j-1,0)-phi(i ,j-1,0)))
1076 + facy*(Real(2.)*(phi(i ,j-1,0)-phi(i ,j ,0))
1077 + (phi(i+1,j-1,0)-phi(i+1,j ,0)))
1078 - fp * (phi(i ,j-1,0)-phi(i ,j ,0)));
1079 }
1080 if (pred(i-1,j)) {
1081 Ax += sig(i-1,j,0)*(facx*(Real(2.)*(phi(i-1,j ,0)-phi(i ,j ,0))
1082 + (phi(i-1,j+1,0)-phi(i ,j+1,0)))
1083 + facy*(Real(2.)*(phi(i ,j+1,0)-phi(i ,j ,0))
1084 + (phi(i-1,j+1,0)-phi(i-1,j ,0)))
1085 + fm * (phi(i ,j+1,0)-phi(i ,j ,0)));
1086 }
1087 if (pred(i,j)) {
1088 Ax += sig(i,j,0)*(facx*(Real(2.)*(phi(i+1,j ,0)-phi(i ,j ,0))
1089 + (phi(i+1,j+1,0)-phi(i ,j+1,0)))
1090 + facy*(Real(2.)*(phi(i ,j+1,0)-phi(i ,j ,0))
1091 + (phi(i+1,j+1,0)-phi(i+1,j ,0)))
1092 - fp * (phi(i ,j+1,0)-phi(i ,j ,0)));
1093 }
1094 return Ax;
1095 }
1096
1097 template <int rr, typename S>
1098 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1099 void mlndlap_Ax_fine_contrib_doit (S const& sig, int i, int j, Box const& ndbx, Box const& ccbx,
1100 Array4<Real> const& f, Array4<Real const> const& res,
1101 Array4<Real const> const& rhs,
1102 Array4<Real const> const& phi,
1103 Array4<int const> const& msk, bool is_rz,
1104 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1105 {
1106 const int ii = rr*i;
1107 const int jj = rr*j;
1108 if (msk(ii,jj,0)) {
1109 Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1110 Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1111
1112 auto is_fine = [&ccbx] (int ix, int iy) -> bool {
1113 return ccbx.contains(ix,iy,0);
1114 };
1115
1116 Real Df = Real(0.0);
1117
1118 const int ilo = amrex::max(ii-rr+1, ndbx.smallEnd(0));
1119 const int ihi = amrex::min(ii+rr-1, ndbx.bigEnd (0));
1120 const int jlo = amrex::max(jj-rr+1, ndbx.smallEnd(1));
1121 const int jhi = amrex::min(jj+rr-1, ndbx.bigEnd (1));
1122
1123 for (int joff = jlo; joff <= jhi; ++joff) {
1124 for (int ioff = ilo; ioff <= ihi; ++ioff) {
1125 Real scale = static_cast<Real>((rr-std::abs(ii-ioff)) *
1126 (rr-std::abs(jj-joff)));
1127 if (ndbx.strictly_contains(ioff,joff,0)) {
1128 Df += scale * (rhs(ioff,joff,0)-res(ioff,joff,0));
1129 } else {
1130 Df += scale * mlndlap_sum_Ax
1131 (is_fine, sig, ioff, joff, facx, facy, phi, is_rz);
1132 }
1133 }}
1134
1135 f(i,j,0) = Df * (Real(1.0)/static_cast<Real>(rr*rr*rr*rr));
1136 } else {
1137 f(i,j,0) = Real(0.0);
1138 }
1139 }
1140
1141template <int rr>
1142AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1143void mlndlap_Ax_fine_contrib (int i, int j, int /*k*/, Box const& ndbx, Box const& ccbx,
1144 Array4<Real> const& f, Array4<Real const> const& res,
1145 Array4<Real const> const& rhs, Array4<Real const> const& phi,
1146 Array4<Real const> const& sig, Array4<int const> const& msk,
1147 bool is_rz,
1148 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1149{
1150 mlndlap_Ax_fine_contrib_doit<rr>
1151 ([&sig] (int ix, int iy, int) -> Real const& { return sig(ix,iy,0); },
1152 i,j,ndbx,ccbx,f,res,rhs,phi,msk,is_rz,dxinv);
1153}
1154
1155template <int rr>
1156AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1157void mlndlap_Ax_fine_contrib_cs (int i, int j, int /*k*/, Box const& ndbx, Box const& ccbx,
1158 Array4<Real> const& f, Array4<Real const> const& res,
1159 Array4<Real const> const& rhs, Array4<Real const> const& phi,
1160 Real const sig, Array4<int const> const& msk,
1161 bool is_rz,
1162 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1163{
1164 mlndlap_Ax_fine_contrib_doit<rr>
1165 ([=] (int, int, int) -> Real { return sig; },
1166 i,j,ndbx,ccbx,f,res,rhs,phi,msk,is_rz,dxinv);
1167}
1168
1169AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1170void mlndlap_res_cf_contrib (int i, int j, int, Array4<Real> const& res,
1171 Array4<Real const> const& phi, Array4<Real const> const& rhs,
1172 Array4<Real const> const& sig, Array4<int const> const& dmsk,
1173 Array4<int const> const& ndmsk, Array4<int const> const& ccmsk,
1174 Array4<Real const> const& fc,
1175 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1176 Box const& ccdom_p, Box const& nddom,
1177 bool is_rz,
1178 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
1179 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
1180 bool neumann_doubling) noexcept
1181{
1182 using namespace nodelap_detail;
1183
1184 if (!dmsk(i,j,0) && ndmsk(i,j,0) == crse_fine_node) {
1185 Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1186 Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1187
1188 Real Ax = mlndlap_sum_Ax([&ccmsk, &ccdom_p] (int ix, int iy) -> bool
1189 {
1190 return ccdom_p.contains(ix,iy,0)
1191 && (ccmsk(ix,iy,0) == crse_cell);
1192 },
1193 [&sig] (int ix, int iy, int) -> Real const&
1194 {
1195 return sig(ix,iy,0);
1196 },
1197 i, j, facx, facy, phi, is_rz);
1198 Ax += fc(i,j,0);
1199 Real const ns = (neumann_doubling) ? neumann_scale(i,j,nddom,bclo,bchi) : Real(1.0);
1200 res(i,j,0) = rhs(i,j,0) - Ax*ns;
1201 }
1202}
1203
1204AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1205void mlndlap_res_cf_contrib_cs (int i, int j, int, Array4<Real> const& res,
1206 Array4<Real const> const& phi, Array4<Real const> const& rhs,
1207 Real const sig, Array4<int const> const& dmsk,
1208 Array4<int const> const& ndmsk, Array4<int const> const& ccmsk,
1209 Array4<Real const> const& fc,
1210 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1211 Box const& ccdom_p, Box const& nddom,
1212 bool is_rz,
1213 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
1214 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi,
1215 bool neumann_doubling) noexcept
1216{
1217 using namespace nodelap_detail;
1218
1219 if (!dmsk(i,j,0) && ndmsk(i,j,0) == crse_fine_node) {
1220 Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1221 Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1222
1223 Real Ax = mlndlap_sum_Ax([&ccmsk, &ccdom_p] (int ix, int iy) -> bool
1224 {
1225 return ccdom_p.contains(ix,iy,0)
1226 && (ccmsk(ix,iy,0) == crse_cell);
1227 },
1228 [=] (int, int, int) -> Real
1229 {
1230 return sig;
1231 },
1232 i, j, facx, facy, phi, is_rz);
1233 Ax += fc(i,j,0);
1234 Real const ns = (neumann_doubling) ? neumann_scale(i,j,nddom,bclo,bchi) : Real(1.0);
1235 res(i,j,0) = rhs(i,j,0) - Ax*ns;
1236 }
1237}
1238
1239//
1240// RAP
1241//
1242
1243AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1244void mlndlap_set_stencil (Box const& bx, Array4<Real> const& sten,
1245 Array4<Real const> const& sigma,
1246 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1247{
1248 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
1249 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
1250 Real fxy = facx + facy;
1251 Real f2xmy = Real(2.0)*facx - facy;
1252 Real fmx2y = Real(2.0)*facy - facx;
1253
1254 amrex::LoopConcurrent(bx, [=] (int i, int j, int k) noexcept
1255 {
1256 sten(i,j,k,1) = f2xmy*(sigma(i,j-1,k)+sigma(i,j,k));
1257 sten(i,j,k,2) = fmx2y*(sigma(i-1,j,k)+sigma(i,j,k));
1258 sten(i,j,k,3) = fxy*sigma(i,j,k);
1259 });
1260}
1261
1262AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1263void mlndlap_set_stencil_s0 (int i, int j, int k, Array4<Real> const& sten) noexcept
1264{
1265 using namespace nodelap_detail;
1266
1267 sten(i,j,k,0) = -(sten(i-1,j ,k,1) + sten(i ,j ,k,1)
1268 + sten(i ,j-1,k,2) + sten(i ,j ,k,2)
1269 + sten(i-1,j-1,k,3) + sten(i ,j-1,k,3)
1270 + sten(i-1,j ,k,3) + sten(i ,j ,k,3));
1271 sten(i,j,k,4) = Real(1.0) / (std::abs(sten(i-1,j ,k,1)) + std::abs(sten(i,j ,k,1))
1272 + std::abs(sten(i ,j-1,k,2)) + std::abs(sten(i,j ,k,2))
1273 + std::abs(sten(i-1,j-1,k,3)) + std::abs(sten(i,j-1,k,3))
1274 + std::abs(sten(i-1,j ,k,3)) + std::abs(sten(i,j ,k,3)) + eps);
1275}
1276
1277AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1278void mlndlap_stencil_rap (int i, int j, int, Array4<Real> const& csten,
1279 Array4<Real const> const& fsten) noexcept
1280{
1281 using namespace nodelap_detail;
1282
1283 constexpr int k = 0;
1284
1285#if 0
1286 auto interp_from_mm_to = [&fsten] (int i_, int j_) -> Real {
1287 Real wxm = std::abs(fsten(i_-1,j_ ,0,1))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_-1,j_ ,0,3))+eps);
1288 Real wym = std::abs(fsten(i_ ,j_-1,0,2))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_ ,j_-1,0,3))+eps);
1289 Real wmm = std::abs(fsten(i_-1,j_-1,0,3)) * (Real(1.) + wxm + wym);
1290 return wmm * fsten(i_,j_,0,4);
1291 };
1292#endif
1293
1294 auto interp_from_mp_to = [&fsten] (int i_, int j_) -> Real {
1295 Real wxm = std::abs(fsten(i_-1,j_ ,0,1))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_-1,j_ ,0,3))+eps);
1296 Real wyp = std::abs(fsten(i_ ,j_ ,0,2))/(std::abs(fsten(i_-1,j_ ,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1297 Real wmp = std::abs(fsten(i_-1,j_ ,0,3)) *(Real(1.) + wxm + wyp);
1298 return wmp * fsten(i_,j_,0,4);
1299 };
1300
1301 auto interp_from_pm_to = [&fsten] (int i_, int j_) -> Real {
1302 Real wxp = std::abs(fsten(i_ ,j_ ,0,1))/(std::abs(fsten(i_ ,j_-1,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1303 Real wym = std::abs(fsten(i_ ,j_-1,0,2))/(std::abs(fsten(i_-1,j_-1,0,3))+std::abs(fsten(i_ ,j_-1,0,3))+eps);
1304 Real wpm = std::abs(fsten(i_ ,j_-1,0,3)) * (Real(1.) + wxp + wym);
1305 return wpm * fsten(i_,j_,0,4);
1306 };
1307
1308 auto interp_from_pp_to = [&fsten] (int i_, int j_) -> Real {
1309 Real wxp = std::abs(fsten(i_ ,j_ ,0,1))/(std::abs(fsten(i_ ,j_-1,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1310 Real wyp = std::abs(fsten(i_ ,j_ ,0,2))/(std::abs(fsten(i_-1,j_ ,0,3))+std::abs(fsten(i_ ,j_ ,0,3))+eps);
1311 Real wpp = std::abs(fsten(i_ ,j_ ,0,3)) * (Real(1.) + wxp + wyp);
1312 return wpp * fsten(i_,j_,0,4);
1313 };
1314
1315 auto interp_from_m0_to = [&fsten] (int i_, int j_) -> Real {
1316 return std::abs(fsten(i_-1,j_,0,1))/(std::abs(fsten(i_-1,j_,0,1))+std::abs(fsten(i_,j_,0,1))+eps);
1317 };
1318
1319 auto interp_from_p0_to = [&fsten] (int i_, int j_) -> Real {
1320 return std::abs(fsten(i_,j_,0,1))/(std::abs(fsten(i_-1,j_,0,1))+std::abs(fsten(i_,j_,0,1))+eps);
1321 };
1322
1323 auto interp_from_0m_to = [&fsten] (int i_, int j_) -> Real {
1324 return std::abs(fsten(i_,j_-1,0,2))/(std::abs(fsten(i_,j_-1,0,2))+std::abs(fsten(i_,j_,0,2))+eps);
1325 };
1326
1327 auto interp_from_0p_to = [&fsten] (int i_, int j_) -> Real {
1328 return std::abs(fsten(i_,j_,0,2))/(std::abs(fsten(i_,j_-1,0,2))+std::abs(fsten(i_,j_,0,2))+eps);
1329 };
1330
1331#if 0
1332 auto Amm = [&fsten] (int i_, int j_) -> Real {
1333 return fsten(i_-1,j_-1,0,3);
1334 };
1335#endif
1336
1337 auto A0m = [&fsten] (int i_, int j_) -> Real {
1338 return fsten(i_,j_-1,0,2);
1339 };
1340
1341 auto Apm = [&fsten] (int i_, int j_) -> Real {
1342 return fsten(i_,j_-1,0,3);
1343 };
1344
1345 auto Am0 = [&fsten] (int i_, int j_) -> Real {
1346 return fsten(i_-1,j_,0,1);
1347 };
1348
1349 auto A00 = [&fsten] (int i_, int j_) -> Real {
1350 return fsten(i_,j_,0,0);
1351 };
1352
1353 auto Ap0 = [&fsten] (int i_, int j_) -> Real {
1354 return fsten(i_,j_,0,1);
1355 };
1356
1357 auto Amp = [&fsten] (int i_, int j_) -> Real {
1358 return fsten(i_-1,j_,0,3);
1359 };
1360
1361 auto A0p = [&fsten] (int i_, int j_) -> Real {
1362 return fsten(i_,j_,0,2);
1363 };
1364
1365 auto App = [&fsten] (int i_, int j_) -> Real {
1366 return fsten(i_,j_,0,3);
1367 };
1368
1369#if 0
1370 auto restrict_from_mm_to = [&fsten] (int ii_, int jj_) -> Real {
1371 Real wxp = std::abs(fsten(ii_-1,jj_-1,0,1))/(std::abs(fsten(ii_-1,jj_-2,0,3))+std::abs(fsten(ii_-1,jj_-1,0,3))+eps);
1372 Real wyp = std::abs(fsten(ii_-1,jj_-1,0,2))/(std::abs(fsten(ii_-2,jj_-1,0,3))+std::abs(fsten(ii_-1,jj_-1,0,3))+eps);
1373 Real wpp = std::abs(fsten(ii_-1,jj_-1,0,3))*(Real(1.)+wxp+wyp);
1374 return wpp * fsten(ii_-1,jj_-1,0,4);
1375 };
1376#endif
1377
1378 auto restrict_from_0m_to = [&fsten] (int ii_, int jj_) -> Real {
1379 return std::abs(fsten(ii_,jj_-1,0,2))/(std::abs(fsten(ii_,jj_-2,0,2))+std::abs(fsten(ii_,jj_-1,0,2))+eps);
1380 };
1381
1382 auto restrict_from_pm_to = [&fsten] (int ii_, int jj_) -> Real {
1383 Real wxm = std::abs(fsten(ii_ ,jj_-1,0,1))/(std::abs(fsten(ii_,jj_-2,0,3))+std::abs(fsten(ii_ ,jj_-1,0,3))+eps);
1384 Real wyp = std::abs(fsten(ii_+1,jj_-1,0,2))/(std::abs(fsten(ii_,jj_-1,0,3))+std::abs(fsten(ii_+1,jj_-1,0,3))+eps);
1385 Real wmp = std::abs(fsten(ii_ ,jj_-1,0,3)) *(Real(1.) + wxm + wyp);
1386 return wmp * fsten(ii_+1,jj_-1,0,4);
1387 };
1388
1389 auto restrict_from_m0_to = [&fsten] (int ii_, int jj_) -> Real {
1390 return std::abs(fsten(ii_-1,jj_,0,1))/(std::abs(fsten(ii_-2,jj_,0,1))+std::abs(fsten(ii_-1,jj_,0,1))+eps);
1391 };
1392
1393 auto restrict_from_p0_to = [&fsten] (int ii_, int jj_) -> Real {
1394 return std::abs(fsten(ii_,jj_,0,1))/(std::abs(fsten(ii_,jj_,0,1))+std::abs(fsten(ii_+1,jj_,0,1))+eps);
1395 };
1396
1397 auto restrict_from_mp_to = [&fsten] (int ii_, int jj_) -> Real {
1398 Real wxp = std::abs(fsten(ii_-1,jj_+1,0,1))/(std::abs(fsten(ii_-1,jj_,0,3))+std::abs(fsten(ii_-1,jj_+1,0,3))+eps);
1399 Real wym = std::abs(fsten(ii_-1,jj_ ,0,2))/(std::abs(fsten(ii_-2,jj_,0,3))+std::abs(fsten(ii_-1,jj_ ,0,3))+eps);
1400 Real wpm = std::abs(fsten(ii_-1,jj_ ,0,3)) * (Real(1.) + wxp + wym);
1401 return wpm * fsten(ii_-1,jj_+1,0,4);
1402 };
1403
1404 auto restrict_from_0p_to = [&fsten] (int ii_, int jj_) -> Real {
1405 return std::abs(fsten(ii_,jj_,0,2))/(std::abs(fsten(ii_,jj_,0,2))+std::abs(fsten(ii_,jj_+1,0,2))+eps);
1406 };
1407
1408 auto restrict_from_pp_to = [&fsten] (int ii_, int jj_) -> Real {
1409 Real wxm = std::abs(fsten(ii_ ,jj_+1,0,1))/(std::abs(fsten(ii_ ,jj_ ,0,3))+std::abs(fsten(ii_ ,jj_+1,0,3))+eps);
1410 Real wym = std::abs(fsten(ii_+1,jj_ ,0,2))/(std::abs(fsten(ii_ ,jj_ ,0,3))+std::abs(fsten(ii_+1,jj_ ,0,3))+eps);
1411 Real wmm = std::abs(fsten(ii_ ,jj_ ,0,3)) * (Real(1.) + wxm + wym);
1412 return wmm * fsten(ii_+1,jj_+1,0,4);
1413 };
1414
1415 int ii = 2*i;
1416 int jj = 2*j;
1417 Array2D<Real,-1,1,-1,1> ap, p;
1418
1419 // csten(i,j,k,1)
1420 p(-1,-1) = interp_from_pp_to(ii+1,jj-1);
1421 p( 0,-1) = interp_from_0p_to(ii+2,jj-1);
1422 p(-1, 0) = interp_from_p0_to(ii+1,jj );
1423 p( 0, 0) = Real(1.);
1424 p(-1, 1) = interp_from_pm_to(ii+1,jj+1);
1425 p( 0, 1) = interp_from_0m_to(ii+2,jj+1);
1426
1427 ap(0,-1) = Ap0(ii,jj-1)*p(-1,-1) + App(ii,jj-1)*p(-1,0);
1428 ap(1,-1) = A00(ii+1,jj-1)*p(-1,-1) + Ap0(ii+1,jj-1)*p(0,-1)
1429 + A0p(ii+1,jj-1)*p(-1,0) + App(ii+1,jj-1)*p(0,0);
1430 ap(0,0) = Apm(ii,jj)*p(-1,-1) + Ap0(ii,jj)*p(-1,0) + App(ii,jj)*p(-1,1);
1431 ap(1,0) = A0m(ii+1,jj)*p(-1,-1) + Apm(ii+1,jj)*p(0,-1)
1432 + A00(ii+1,jj)*p(-1,0) + Ap0(ii+1,jj)*p(0,0)
1433 + A0p(ii+1,jj)*p(-1,1) + App(ii+1,jj)*p(0,1);
1434 ap(0,1) = Apm(ii,jj+1)*p(-1,0) + Ap0(ii,jj+1)*p(-1,1);
1435 ap(1,1) = A0m(ii+1,jj+1)*p(-1,0) + Apm(ii+1,jj+1)*p(0,0)
1436 + A00(ii+1,jj+1)*p(-1,1) + Ap0(ii+1,jj+1)*p(0,1);
1437
1438 csten(i,j,k,1) = Real(0.25)*(restrict_from_0m_to(ii,jj)*ap(0,-1)
1439 + restrict_from_pm_to(ii,jj)*ap(1,-1)
1440 + ap(0,0)
1441 + restrict_from_p0_to(ii,jj)*ap(1,0)
1442 + restrict_from_0p_to(ii,jj)*ap(0,1)
1443 + restrict_from_pp_to(ii,jj)*ap(1,1));
1444
1445 // csten(i,j,k,2)
1446 p(-1,-1) = interp_from_pp_to(ii-1,jj+1);
1447 p( 0,-1) = interp_from_0p_to(ii ,jj+1);
1448 p( 1,-1) = interp_from_mp_to(ii+1,jj+1);
1449 p(-1, 0) = interp_from_p0_to(ii-1,jj+2);
1450 p( 0, 0) = Real(1.);
1451 p( 1, 0) = interp_from_m0_to(ii+1,jj+2);
1452
1453 ap(-1,0) = A0p(ii-1,jj)*p(-1,-1) + App(ii-1,jj)*p(0,-1);
1454 ap(0,0) = Amp(ii,jj)*p(-1,-1) + A0p(ii,jj)*p(0,-1) + App(ii,jj)*p(1,-1);
1455 ap(1,0) = Amp(ii+1,jj)*p(0,-1) + A0p(ii+1,jj)*p(1,-1);
1456 ap(-1,1) = A00(ii-1,jj+1)*p(-1,-1) + Ap0(ii-1,jj+1)*p(0,-1)
1457 + A0p(ii-1,jj+1)*p(-1,0) + App(ii-1,jj+1)*p(0,0);
1458 ap(0,1) = Am0(ii,jj+1)*p(-1,-1) + A00(ii,jj+1)*p(0,-1) + Ap0(ii,jj+1)*p(1,-1)
1459 + Amp(ii,jj+1)*p(-1,0) + A0p(ii,jj+1)*p(0,0) + App(ii,jj+1)*p(1,0);
1460 ap(1,1) = Am0(ii+1,jj+1)*p(0,-1) + A00(ii+1,jj+1)*p(1,-1)
1461 + Amp(ii+1,jj+1)*p(0,0) + A0p(ii+1,jj+1)*p(1,0);
1462
1463 csten(i,j,k,2) = Real(0.25)*(restrict_from_m0_to(ii,jj)*ap(-1,0)
1464 + ap(0,0)
1465 + restrict_from_p0_to(ii,jj)*ap(1,0)
1466 + restrict_from_mp_to(ii,jj)*ap(-1,1)
1467 + restrict_from_0p_to(ii,jj)*ap(0,1)
1468 + restrict_from_pp_to(ii,jj)*ap(1,1));
1469
1470 // csten(i,j,k,3)
1471 p(-1,-1) = interp_from_pp_to(ii+1,jj+1);
1472 p( 0,-1) = interp_from_0p_to(ii+2,jj+1);
1473 p(-1, 0) = interp_from_p0_to(ii+1,jj+2);
1474 p( 0, 0) = Real(1.);
1475
1476 ap(0,0) = App(ii,jj)*p(-1,-1);
1477 ap(1,0) = A0p(ii+1,jj)*p(-1,-1) + App(ii+1,jj)*p(0,-1);
1478 ap(0,1) = Ap0(ii,jj+1)*p(-1,-1) + App(ii,jj+1)*p(-1,0);
1479 ap(1,1) = A00(ii+1,jj+1)*p(-1,-1) + Ap0(ii+1,jj+1)*p(0,-1)
1480 + A0p(ii+1,jj+1)*p(-1,0) + App(ii+1,jj+1)*p(0,0);
1481
1482 Real cross1 = Real(0.25)*(ap(0,0)
1483 + restrict_from_p0_to(ii,jj)*ap(1,0)
1484 + restrict_from_0p_to(ii,jj)*ap(0,1)
1485 + restrict_from_pp_to(ii,jj)*ap(1,1));
1486
1487 p(0,-1) = interp_from_0p_to(ii,jj+1);
1488 p(1,-1) = interp_from_mp_to(ii+1,jj+1);
1489 p(0, 0) = Real(1.);
1490 p(1, 0) = interp_from_m0_to(ii+1,jj+2);
1491
1492 ap(-1,0) = Amp(ii+1,jj)*p(0,-1) + A0p(ii+1,jj)*p(1,-1);
1493 ap( 0,0) = Amp(ii+2,jj)*p(1,-1);
1494 ap(-1,1) = Am0(ii+1,jj+1)*p(0,-1) + A00(ii+1,jj+1)*p(1,-1) + Amp(ii+1,jj+1)*p(0,0)
1495 + A0p(ii+1,jj+1)*p(1,0);
1496 ap( 0,1) = Am0(ii+2,jj+1)*p(1,-1) + Amp(ii+2,jj+1)*p(1,0);
1497
1498 Real cross2 = Real(0.25)*(ap(0,0)
1499 + restrict_from_m0_to(ii+2,jj)*ap(-1,0)
1500 + restrict_from_mp_to(ii+2,jj)*ap(-1,1)
1501 + restrict_from_0p_to(ii+2,jj)*ap( 0,1));
1502
1503 csten(i,j,k,3) = Real(0.5)*(cross1+cross2);
1504}
1505
1506AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1507Real mlndlap_adotx_sten_doit (int i, int j, int k, Array4<Real const> const& x,
1508 Array4<Real const> const& sten) noexcept
1509{
1510 return x(i-1,j-1,k)*sten(i-1,j-1,k,3)
1511 + x(i ,j-1,k)*sten(i ,j-1,k,2)
1512 + x(i+1,j-1,k)*sten(i ,j-1,k,3)
1513 + x(i-1,j ,k)*sten(i-1,j ,k,1)
1514 + x(i ,j ,k)*sten(i ,j ,k,0)
1515 + x(i+1,j ,k)*sten(i ,j ,k,1)
1516 + x(i-1,j+1,k)*sten(i-1,j ,k,3)
1517 + x(i ,j+1,k)*sten(i ,j ,k,2)
1518 + x(i+1,j+1,k)*sten(i ,j ,k,3);
1519}
1520
1521AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1522Real mlndlap_adotx_sten (int i, int j, int k, Array4<Real const> const& x,
1523 Array4<Real const> const& sten, Array4<int const> const& msk) noexcept
1524{
1525 if (msk(i,j,k)) {
1526 return Real(0.0);
1527 } else {
1528 return mlndlap_adotx_sten_doit(i,j,k,x,sten);
1529 }
1530}
1531
1532AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1533void mlndlap_gauss_seidel_sten (int i, int j, int k, Array4<Real> const& sol,
1534 Array4<Real const> const& rhs,
1535 Array4<Real const> const& sten,
1536 Array4<int const> const& msk) noexcept
1537{
1538 if (msk(i,j,k)) {
1539 sol(i,j,k) = Real(0.0);
1540 } else if (sten(i,j,k,0) != Real(0.0)) {
1541 Real Ax = mlndlap_adotx_sten_doit(i,j,k,sol,sten);
1542 sol(i,j,k) += (rhs(i,j,k) - Ax) / sten(i,j,k,0);
1543 }
1544}
1545
1546inline
1547void mlndlap_gauss_seidel_sten (Box const& bx, Array4<Real> const& sol,
1548 Array4<Real const> const& rhs,
1549 Array4<Real const> const& sten,
1550 Array4<int const> const& msk) noexcept
1551{
1552 AMREX_LOOP_3D(bx, i, j, k,
1553 {
1554 mlndlap_gauss_seidel_sten(i,j,k,sol,rhs,sten,msk);
1555 });
1556}
1557
1558AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1559void mlndlap_interpadd_rap (int i, int j, int, Array4<Real> const& fine,
1560 Array4<Real const> const& crse, Array4<Real const> const& sten,
1561 Array4<int const> const& msk) noexcept
1562{
1563 using namespace nodelap_detail;
1564
1565 if (!msk(i,j,0) && sten(i,j,0,0) != Real(0.0)) {
1566 int ic = amrex::coarsen(i,2);
1567 int jc = amrex::coarsen(j,2);
1568 bool ieven = ic*2 == i;
1569 bool jeven = jc*2 == j;
1570 Real fv;
1571 if (ieven && jeven) {
1572 fv = crse(ic,jc,0);
1573 } else if (ieven) {
1574 Real wym = std::abs(sten(i,j-1,0,2));
1575 Real wyp = std::abs(sten(i,j ,0,2));
1576 fv = (wym*crse(ic,jc,0) + wyp*crse(ic,jc+1,0)) / (wym+wyp+eps);
1577 } else if (jeven) {
1578 Real wxm = std::abs(sten(i-1,j,0,1));
1579 Real wxp = std::abs(sten(i ,j,0,1));
1580 fv = (wxm*crse(ic,jc,0) + wxp*crse(ic+1,jc,0)) / (wxm+wxp+eps);
1581 } else {
1582 Real wxm = std::abs(sten(i-1,j ,0,1)) /
1583 (std::abs(sten(i-1,j-1,0,3))+std::abs(sten(i-1,j ,0,3))+eps);
1584 Real wxp = std::abs(sten(i ,j ,0,1)) /
1585 (std::abs(sten(i ,j-1,0,3))+std::abs(sten(i ,j ,0,3))+eps);
1586 Real wym = std::abs(sten(i ,j-1,0,2)) /
1587 (std::abs(sten(i-1,j-1,0,3))+std::abs(sten(i ,j-1,0,3))+eps);
1588 Real wyp = std::abs(sten(i ,j ,0,2)) /
1589 (std::abs(sten(i-1,j ,0,3))+std::abs(sten(i ,j ,0,3))+eps);
1590 Real wmm = std::abs(sten(i-1,j-1,0,3)) * (Real(1.0) + wxm + wym);
1591 Real wpm = std::abs(sten(i,j-1,0,3)) * (Real(1.0) + wxp + wym);
1592 Real wmp = std::abs(sten(i-1,j,0,3)) *(Real(1.0) + wxm + wyp);
1593 Real wpp = std::abs(sten(i,j,0,3)) * (Real(1.0) + wxp + wyp);
1594 fv = (wmm*crse(ic,jc,0) + wpm*crse(ic+1,jc,0)
1595 + wmp*crse(ic,jc+1,0) + wpp*crse(ic+1,jc+1,0))
1596 / (wmm+wpm+wmp+wpp+eps);
1597 }
1598
1599 fine(i,j,0) += fv;
1600 }
1601}
1602
1603AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1604void mlndlap_restriction_rap (int i, int j, int /*k*/, Array4<Real> const& crse,
1605 Array4<Real const> const& fine, Array4<Real const> const& sten,
1606 Array4<int const> const& msk) noexcept
1607{
1608 using namespace nodelap_detail;
1609
1610 int ii = i*2;
1611 int jj = j*2;
1612 if (msk(ii,jj,0)) {
1613 crse(i,j,0) = Real(0.0);
1614 } else {
1615
1616 Real cv = fine(ii,jj,0)
1617 + fine(ii-1,jj ,0)*std::abs(sten(ii-1,jj ,0,1))
1618 / (std::abs(sten(ii-2,jj ,0,1))
1619 +std::abs(sten(ii-1,jj ,0,1))+eps)
1620 + fine(ii+1,jj ,0)*std::abs(sten(ii ,jj ,0,1))
1621 / (std::abs(sten(ii ,jj ,0,1))
1622 +std::abs(sten(ii+1,jj ,0,1))+eps)
1623 + fine(ii ,jj-1,0)*std::abs(sten(ii ,jj-1,0,2))
1624 / (std::abs(sten(ii ,jj-2,0,2))
1625 +std::abs(sten(ii ,jj-1,0,2))+eps)
1626 + fine(ii ,jj+1,0)*std::abs(sten(ii ,jj ,0,2))
1627 / (std::abs(sten(ii ,jj ,0,2))
1628 +std::abs(sten(ii ,jj+1,0,2))+eps);
1629
1630 Real wxp = std::abs(sten(ii-1,jj-1,0,1))
1631 / (std::abs(sten(ii-1,jj-2,0,3))
1632 +std::abs(sten(ii-1,jj-1,0,3))+eps);
1633 Real wyp = std::abs(sten(ii-1,jj-1,0,2))
1634 / (std::abs(sten(ii-2,jj-1,0,3))
1635 +std::abs(sten(ii-1,jj-1,0,3))+eps);
1636 Real wpp = std::abs(sten(ii-1,jj-1,0,3))*(Real(1.0) + wxp + wyp);
1637 cv += wpp*sten(ii-1,jj-1,0,4)*fine(ii-1,jj-1,0);
1638
1639 Real wxm = std::abs(sten(ii ,jj-1,0,1))
1640 / (std::abs(sten(ii ,jj-2,0,3))
1641 +std::abs(sten(ii ,jj-1,0,3))+eps);
1642 wyp = std::abs(sten(ii+1,jj-1,0,2))
1643 / (std::abs(sten(ii ,jj-1,0,3))
1644 +std::abs(sten(ii+1,jj-1,0,3))+eps);
1645 Real wmp = std::abs(sten(ii ,jj-1,0,3))*(Real(1.0) + wxm + wyp);
1646 cv += wmp*sten(ii+1,jj-1,0,4)*fine(ii+1,jj-1,0);
1647
1648 wxp = std::abs(sten(ii-1,jj+1,0,1))
1649 / (std::abs(sten(ii-1,jj ,0,3))
1650 +std::abs(sten(ii-1,jj+1,0,3))+eps);
1651 Real wym = std::abs(sten(ii-1,jj ,0,2))
1652 / (std::abs(sten(ii-2,jj ,0,3))
1653 +std::abs(sten(ii-1,jj ,0,3))+eps);
1654 Real wpm = std::abs(sten(ii-1,jj ,0,3)) * (Real(1.0) + wxp + wym);
1655 cv += wpm*sten(ii-1,jj+1,0,4)*fine(ii-1,jj+1,0);
1656
1657 wxm = std::abs(sten(ii ,jj+1,0,1))
1658 / (std::abs(sten(ii ,jj ,0,3))
1659 +std::abs(sten(ii ,jj+1,0,3))+eps);
1660 wym = std::abs(sten(ii+1,jj ,0,2))
1661 / (std::abs(sten(ii ,jj ,0,3))
1662 +std::abs(sten(ii+1,jj ,0,3))+eps);
1663 Real wmm = std::abs(sten(ii ,jj ,0,3)) * (Real(1.0) + wxm + wym);
1664 cv += wmm*sten(ii+1,jj+1,0,4)*fine(ii+1,jj+1,0);
1665
1666 crse(i,j,0) = cv * Real(0.25);
1667 }
1668}
1669
1670#ifdef AMREX_USE_EB
1671
1673namespace nodelap_detail {
1674 constexpr int i_S_x = 0;
1675 constexpr int i_S_y = 1;
1676 constexpr int i_S_x2 = 2;
1677 constexpr int i_S_y2 = 3;
1678 constexpr int i_S_xy = 4;
1679 constexpr int n_Sintg = 5;
1680
1681 constexpr int i_B_x = 0;
1682 constexpr int i_B_y = 1;
1683 constexpr int i_B_xy = 2;
1684 constexpr int n_Bintg = 3;
1685}
1687
1688AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1689void mlndlap_set_connection (int i, int j, int, Array4<Real> const& conn,
1690 Array4<Real const> const& intg, Array4<Real const> const& vol,
1691 Array4<EBCellFlag const> const& flag) noexcept
1692{
1693 using namespace nodelap_detail;
1694
1695 if (flag(i,j,0).isCovered()) {
1696 for (int n = 0; n < 6; ++n) { conn(i,j,0,n) = Real(0.); }
1697 } else if (flag(i,j,0).isRegular() || vol(i,j,0) >= almostone) {
1698 for (int n = 0; n < 6; ++n) { conn(i,j,0,n) = Real(1.); }
1699 } else {
1700 // Note that these are normalized so that they equal 1 in the case of a regular cell
1701
1702 conn(i,j,0,0) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_y2) - intg(i,j,0,i_S_y));
1703 conn(i,j,0,1) = Real(6.)*(Real(0.25)*vol(i,j,0) - intg(i,j,0,i_S_y2));
1704 conn(i,j,0,2) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_y2) + intg(i,j,0,i_S_y));
1705
1706 conn(i,j,0,3) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_x2) - intg(i,j,0,i_S_x));
1707 conn(i,j,0,4) = Real(6.)*(Real(0.25)*vol(i,j,0) - intg(i,j,0,i_S_x2));
1708 conn(i,j,0,5) = Real(3.)*(Real(0.25)*vol(i,j,0) + intg(i,j,0,i_S_x2) + intg(i,j,0,i_S_x));
1709 }
1710}
1711
1712AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1713void mlndlap_set_stencil_eb (int i, int j, int, Array4<Real> const& sten,
1714 Array4<Real const> const& sig, Array4<Real const> const& conn,
1715 GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1716{
1717 Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
1718 Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
1719
1720 sten(i,j,0,1) = Real(2.)*facx*(sig(i,j-1,0)*conn(i,j-1,0,2)+sig(i,j,0)*conn(i,j,0,0))
1721 -facy*(sig(i,j-1,0)*conn(i,j-1,0,4)+sig(i,j,0)*conn(i,j,0,4));
1722 sten(i,j,0,2) = Real(2.)*facy*(sig(i-1,j,0)*conn(i-1,j,0,5)+sig(i,j,0)*conn(i,j,0,3))
1723 -facx*(sig(i-1,j,0)*conn(i-1,j,0,1)+sig(i,j,0)*conn(i,j,0,1));
1724 sten(i,j,0,3) = (facx*conn(i,j,0,1)+facy*conn(i,j,0,4))*sig(i,j,0);
1725}
1726
1727
1728AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1729void mlndlap_divu_eb (int i, int j, int, Array4<Real> const& rhs, Array4<Real const> const& vel,
1730 Array4<Real const> const& vfrac, Array4<Real const> const& intg,
1731 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1732 Box const& nodal_domain,
1733 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bclo,
1734 GpuArray<LinOpBCType, AMREX_SPACEDIM> const& bchi) noexcept
1735{
1736 Real facx = Real(0.5)*dxinv[0];
1737 Real facy = Real(0.5)*dxinv[1];
1738
1739 const auto domlo = amrex::lbound(nodal_domain);
1740 const auto domhi = amrex::ubound(nodal_domain);
1741
1742 if (!msk(i,j,0)) {
1743
1744 Real zero_ilo = Real(1.0);
1745 Real zero_ihi = Real(1.0);
1746 Real zero_jlo = Real(1.0);
1747 Real zero_jhi = Real(1.0);
1748
1749 // The nodal divergence operator should not see the tangential velocity
1750 // at an inflow face
1751 if ((bclo[0] == LinOpBCType::Neumann || bclo[0] == LinOpBCType::inflow)
1752 && i == domlo.x)
1753 {
1754 zero_ilo = Real(0.0);
1755 }
1756 if ((bchi[0] == LinOpBCType::Neumann || bchi[0] == LinOpBCType::inflow)
1757 && i == domhi.x)
1758 {
1759 zero_ihi = Real(0.0);
1760 }
1761 if ((bclo[1] == LinOpBCType::Neumann || bclo[1] == LinOpBCType::inflow)
1762 && j == domlo.y)
1763 {
1764 zero_jlo = Real(0.0);
1765 }
1766 if ((bchi[1] == LinOpBCType::Neumann || bchi[1] == LinOpBCType::inflow)
1767 && j == domhi.y)
1768 {
1769 zero_jhi = Real(0.0);
1770 }
1771
1772 rhs(i,j,0) = facx*(-vel(i-1,j-1,0,0)*(vfrac(i-1,j-1,0)+Real(2.)*intg(i-1,j-1,0,1))*zero_jlo
1773 +vel(i ,j-1,0,0)*(vfrac(i ,j-1,0)+Real(2.)*intg(i ,j-1,0,1))*zero_jlo
1774 -vel(i-1,j ,0,0)*(vfrac(i-1,j ,0)-Real(2.)*intg(i-1,j ,0,1))*zero_jhi
1775 +vel(i ,j ,0,0)*(vfrac(i ,j ,0)-Real(2.)*intg(i ,j ,0,1))*zero_jhi)
1776 + facy*(-vel(i-1,j-1,0,1)*(vfrac(i-1,j-1,0)+Real(2.)*intg(i-1,j-1,0,0))*zero_ilo
1777 -vel(i ,j-1,0,1)*(vfrac(i ,j-1,0)-Real(2.)*intg(i ,j-1,0,0))*zero_ihi
1778 +vel(i-1,j ,0,1)*(vfrac(i-1,j ,0)+Real(2.)*intg(i-1,j ,0,0))*zero_ilo
1779 +vel(i ,j ,0,1)*(vfrac(i ,j ,0)-Real(2.)*intg(i ,j ,0,0))*zero_ihi);
1780 } else {
1781 rhs(i,j,0) = Real(0.);
1782 }
1783}
1784
1785AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1786void add_eb_flow_contribution (int i, int j, int, Array4<Real> const& rhs,
1787 Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
1788 Array4<Real const> const& bareaarr,
1789 Array4<Real const> const& sintg,
1790 Array4<Real const> const& eb_vel_dot_n) noexcept
1791{
1792 using namespace nodelap_detail;
1793
1794 Real fac_eb = 0.25 *dxinv[0];
1795
1796 if (!msk(i,j,0)) {
1797 rhs(i,j,0) += fac_eb*(
1798 eb_vel_dot_n(i-1,j-1,0)*( bareaarr(i-1,j-1,0)
1799 +Real(2.)*sintg(i-1,j-1,0,i_B_x )
1800 +Real(2.)*sintg(i-1,j-1,0,i_B_y )
1801 +Real(4.)*sintg(i-1,j-1,0,i_B_xy))
1802 +eb_vel_dot_n(i ,j-1,0)*( bareaarr(i ,j-1,0)
1803 -Real(2.)*sintg(i ,j-1,0,i_B_x )
1804 +Real(2.)*sintg(i ,j-1,0,i_B_y )
1805 -Real(4.)*sintg(i ,j-1,0,i_B_xy))
1806 +eb_vel_dot_n(i-1,j ,0)*( bareaarr(i-1,j ,0)
1807 +Real(2.)*sintg(i-1,j ,0,i_B_x )
1808 -Real(2.)*sintg(i-1,j ,0,i_B_y )
1809 -Real(4.)*sintg(i-1,j ,0,i_B_xy))
1810 +eb_vel_dot_n(i ,j ,0)*( bareaarr(i ,j ,0)
1811 -Real(2.)*sintg(i ,j ,0,i_B_x )
1812 -Real(2.)*sintg(i ,j ,0,i_B_y )
1813 +Real(4.)*sintg(i ,j ,0,i_B_xy)));
1814 }
1815}
1816
1817AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1818void mlndlap_mknewu_eb (int i, int j, int, Array4<Real> const& u, Array4<Real const> const& p,
1819 Array4<Real const> const& sig, Array4<Real const> const& vfrac,
1820 Array4<Real const> const& intg, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1821{
1822 Real facx = Real(0.5)*dxinv[0];
1823 Real facy = Real(0.5)*dxinv[1];
1824 if (vfrac(i,j,0) == Real(0.)) {
1825 u(i,j,0,0) = u(i,j,0,1) = Real(0.);
1826 } else {
1827 Real dpdx = facx*(-p(i,j,0)+p(i+1,j,0)-p(i,j+1,0)+p(i+1,j+1,0));
1828 Real dpdy = facy*(-p(i,j,0)-p(i+1,j,0)+p(i,j+1,0)+p(i+1,j+1,0));
1829 Real dpp = (p(i,j,0)+p(i+1,j+1,0)-p(i+1,j,0)-p(i,j+1,0))/vfrac(i,j,0);
1830 u(i,j,0,0) -= sig(i,j,0)*(dpdx + dxinv[0]*intg(i,j,0,1)*dpp);
1831 u(i,j,0,1) -= sig(i,j,0)*(dpdy + dxinv[1]*intg(i,j,0,0)*dpp);
1832 }
1833}
1834
1835AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1836void mlndlap_mknewu_eb_c (int i, int j, int, Array4<Real> const& u, Array4<Real const> const& p,
1837 Real sig, Array4<Real const> const& vfrac,
1838 Array4<Real const> const& intg, GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
1839{
1840 Real facx = Real(0.5)*dxinv[0];
1841 Real facy = Real(0.5)*dxinv[1];
1842 if (vfrac(i,j,0) == Real(0.)) {
1843 u(i,j,0,0) = u(i,j,0,1) = Real(0.);
1844 } else {
1845 Real dpdx = facx*(-p(i,j,0)+p(i+1,j,0)-p(i,j+1,0)+p(i+1,j+1,0));
1846 Real dpdy = facy*(-p(i,j,0)-p(i+1,j,0)+p(i,j+1,0)+p(i+1,j+1,0));
1847 Real dpp = (p(i,j,0)+p(i+1,j+1,0)-p(i+1,j,0)-p(i,j+1,0))/vfrac(i,j,0);
1848 u(i,j,0,0) -= sig*(dpdx + dxinv[0]*intg(i,j,0,1)*dpp);
1849 u(i,j,0,1) -= sig*(dpdy + dxinv[1]*intg(i,j,0,0)*dpp);
1850 }
1851}
1852
1853AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1854Real mlndlap_rhcc_eb (int i, int j, int, Array4<Real const> const& rhcc,
1855 Array4<Real const> const& vfrac, Array4<Real const> const& intg,
1856 Array4<int const> const& msk) noexcept
1857{
1858 using namespace nodelap_detail;
1859
1860 if (!msk(i,j,0)) {
1861 return
1862 rhcc(i ,j ,0)*(Real(0.25)*vfrac(i ,j ,0)-intg(i ,j ,0,i_S_x)-intg(i ,j ,0,i_S_y)+intg(i ,j ,0,i_S_xy)) +
1863 rhcc(i-1,j ,0)*(Real(0.25)*vfrac(i-1,j ,0)+intg(i-1,j ,0,i_S_x)-intg(i-1,j ,0,i_S_y)-intg(i-1,j ,0,i_S_xy)) +
1864 rhcc(i-1,j-1,0)*(Real(0.25)*vfrac(i-1,j-1,0)+intg(i-1,j-1,0,i_S_x)+intg(i-1,j-1,0,i_S_y)+intg(i-1,j-1,0,i_S_xy)) +
1865 rhcc(i ,j-1,0)*(Real(0.25)*vfrac(i ,j-1,0)-intg(i ,j-1,0,i_S_x)+intg(i ,j-1,0,i_S_y)-intg(i ,j-1,0,i_S_xy));
1866 } else {
1867 return Real(0.);
1868 }
1869}
1870
1871AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1872void mlndlap_set_integral (int i, int j, int, Array4<Real> const& intg) noexcept
1873{
1874 using namespace nodelap_detail;
1875
1876 intg(i,j,0,i_S_x ) = Real(0.);
1877 intg(i,j,0,i_S_y ) = Real(0.);
1878 intg(i,j,0,i_S_x2) = Real(1./12.);
1879 intg(i,j,0,i_S_y2) = Real(1./12.);
1880 intg(i,j,0,i_S_xy) = Real(0.);
1881}
1882
1883AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1884void mlndlap_set_surface_integral (int i, int j, int, Array4<Real> const& sintg) noexcept
1885{
1886 using namespace nodelap_detail;
1887
1888 sintg(i,j,0,i_B_x ) = Real(0.);
1889 sintg(i,j,0,i_B_y ) = Real(0.);
1890 sintg(i,j,0,i_B_xy) = Real(0.);
1891}
1892
1893AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1894void mlndlap_set_integral_eb (int i, int j, int, Array4<Real> const& intg,
1895 Array4<EBCellFlag const> const& flag, Array4<Real const> const& vol,
1896 Array4<Real const> const& ax, Array4<Real const> const& ay,
1897 Array4<Real const> const& bcen) noexcept
1898{
1899 using namespace nodelap_detail;
1900
1901 if (flag(i,j,0).isCovered()) {
1902 intg(i,j,0,i_S_x ) = Real(0.);
1903 intg(i,j,0,i_S_y ) = Real(0.);
1904 intg(i,j,0,i_S_x2) = Real(0.);
1905 intg(i,j,0,i_S_y2) = Real(0.);
1906 intg(i,j,0,i_S_xy) = Real(0.);
1907 } else if (flag(i,j,0).isRegular() || vol(i,j,0) >= almostone) {
1908 intg(i,j,0,i_S_x ) = Real(0.);
1909 intg(i,j,0,i_S_y ) = Real(0.);
1910 intg(i,j,0,i_S_x2) = Real(1./12.);
1911 intg(i,j,0,i_S_y2) = Real(1./12.);
1912 intg(i,j,0,i_S_xy) = Real(0.);
1913 } else {
1914 Real axm = ax(i,j,0);
1915 Real axp = ax(i+1,j,0);
1916 Real aym = ay(i,j,0);
1917 Real ayp = ay(i,j+1,0);
1918
1919 Real apnorm = std::sqrt((axm-axp)*(axm-axp) + (aym-ayp)*(aym-ayp));
1920 if (apnorm == Real(0.)) {
1921 amrex::Abort("amrex_mlndlap_set_integral: we are in trouble");
1922 }
1923
1924 Real apnorminv = Real(1.)/apnorm;
1925 Real anrmx = (axm-axp) * apnorminv; // pointing to the wall
1926 Real anrmy = (aym-ayp) * apnorminv;
1927
1928 Real bcx = bcen(i,j,0,0);
1929 Real bcy = bcen(i,j,0,1);
1930
1931 Real Sx, Sy, Sx2, Sy2, Sxy;
1932 if (std::abs(anrmx) <= almostzero) {
1933 Sx = Real(0.);
1934 Sx2 = Real(1./24.)*(axm+axp);
1935 Sxy = Real(0.);
1936 } else if (std::abs(anrmy) <= almostzero) {
1937 Sx = Real(1./8.) *(axp-axm) + anrmx*Real(0.5)*(bcx*bcx);
1938 Sx2 = Real(1./24.)*(axp+axm) + anrmx*Real(1./3.)*(bcx*bcx*bcx);
1939 Sxy = Real(0.);
1940 } else {
1941 Real xmin, xmax;
1942 if (anrmx > Real(0.)) {
1943 xmin = Real(-0.5) + amrex::min(aym,ayp);
1944 xmax = Real(-0.5) + amrex::max(aym,ayp);
1945 } else {
1946 xmin = Real(0.5) - amrex::max(aym,ayp);
1947 xmax = Real(0.5) - amrex::min(aym,ayp);
1948 }
1949 Real xmin3 = xmin*xmin*xmin;
1950 Real xmin4 = xmin3*xmin;
1951 Real xmax3 = xmax*xmax*xmax;
1952 Real xmax4 = xmax3*xmax;
1953 Sx = Real(1./8.) *(axp-axm) + (anrmx/std::abs(anrmy))*Real(1./6.) *(xmax3-xmin3);
1954 Sx2 = Real(1./24.)*(axp+axm) + (anrmx/std::abs(anrmy))*Real(1./12.)*(xmax4-xmin4);
1955
1956 Real kk = -anrmx/anrmy;
1957 Real bb = bcy-kk*bcx;
1958 Sxy = Real(1./8.)*kk*kk*(xmax4-xmin4) + Real(1./3.)*kk*bb*(xmax3-xmin3)
1959 + (Real(0.25)*bb*bb-Real(1./16.))*(xmax*xmax-xmin*xmin);
1960 Sxy = std::copysign(Sxy, anrmy);
1961 }
1962
1963 if (std::abs(anrmy) <= almostzero) {
1964 Sy = Real(0.);
1965 Sy2 = Real(1./24.)*(aym+ayp);
1966 } else if (std::abs(anrmx) <= almostzero) {
1967 Sy = Real(1./8.) *(ayp-aym) + anrmy*Real(0.5)*(bcy*bcy);
1968 Sy2 = Real(1./24.)*(ayp+aym) + anrmy*Real(1./3.)*(bcy*bcy*bcy);
1969 } else {
1970 Real ymin, ymax;
1971 if (anrmy > Real(0.)) {
1972 ymin = Real(-0.5) + amrex::min(axm,axp);
1973 ymax = Real(-0.5) + amrex::max(axm,axp);
1974 } else {
1975 ymin = Real(0.5) - amrex::max(axm,axp);
1976 ymax = Real(0.5) - amrex::min(axm,axp);
1977 }
1978 Real ymin3 = ymin*ymin*ymin;
1979 Real ymin4 = ymin3*ymin;
1980 Real ymax3 = ymax*ymax*ymax;
1981 Real ymax4 = ymax3*ymax;
1982 Sy = Real(1./8.) *(ayp-aym) + (anrmy/std::abs(anrmx))*Real(1./6.) *(ymax3-ymin3);
1983 Sy2 = Real(1./24.)*(ayp+aym) + (anrmy/std::abs(anrmx))*Real(1./12.)*(ymax4-ymin4);
1984 }
1985
1986 intg(i,j,0,i_S_x ) = Sx;
1987 intg(i,j,0,i_S_y ) = Sy;
1988 intg(i,j,0,i_S_x2) = Sx2;
1989 intg(i,j,0,i_S_y2) = Sy2;
1990 intg(i,j,0,i_S_xy) = Sxy;
1991 }
1992}
1993
1994AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1995void mlndlap_set_surface_integral_eb (int i, int j, int, Array4<Real> const& sintg,
1996 Array4<EBCellFlag const> const& flag,
1997 Array4<Real const> const& bcen,
1998 Array4<Real const> const& barea,
1999 Array4<Real const> const& bnorm) noexcept
2000{
2001 using namespace nodelap_detail;
2002
2003 if (flag(i,j,0).isCovered() || flag(i,j,0).isRegular()) {
2004 sintg(i,j,0,i_B_x ) = Real(0.);
2005 sintg(i,j,0,i_B_y ) = Real(0.);
2006 sintg(i,j,0,i_B_xy) = Real(0.);
2007 } else {
2008 Real bcx = bcen(i,j,0,0);
2009 Real bcy = bcen(i,j,0,1);
2010
2011 Real btanx = bnorm(i,j,0,1);
2012 Real btany = -bnorm(i,j,0,0);
2013
2014 Real x0 = bcx - Real(0.5)*barea(i,j,0)*btanx;
2015 Real x1 = bcx + Real(0.5)*barea(i,j,0)*btanx;
2016
2017 Real y0 = bcy - Real(0.5)*barea(i,j,0)*btany;
2018 Real y1 = bcy + Real(0.5)*barea(i,j,0)*btany;
2019
2020 Real Bx = barea(i,j,0)*Real(0.5)*(x1 + x0);
2021 Real By = barea(i,j,0)*Real(0.5)*(y1 + y0);
2022 Real Bxy = barea(i,j,0)*(x0*y0 + (x0*(y1 - y0) + y0*(x1 - x0))/Real(2.) + (x1 - x0)*(y1 - y0)/Real(3.));
2023
2024 sintg(i,j,0,i_B_x ) = Bx;
2025 sintg(i,j,0,i_B_y ) = By;
2026 sintg(i,j,0,i_B_xy) = Bxy;
2027 }
2028}
2029
2030#endif
2031
2032#if defined(AMREX_USE_HYPRE)
2033
2034template <typename HypreInt, typename AtomicInt>
2035void mlndlap_fillijmat_sten_cpu (Box const& ndbx,
2036 Array4<AtomicInt const> const& gid,
2037 Array4<int const> const& lid,
2038 HypreInt* ncols, HypreInt* cols,
2039 Real* mat, // NOLINT(readability-non-const-parameter)
2040 Array4<Real const> const& sten) noexcept
2041{
2042 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2043 HypreInt nelems = 0;
2044 amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2045 {
2046 if (lid(i,j,k) >= 0)
2047 {
2048 cols[nelems] = gid(i,j,k);
2049 mat[nelems] = sten(i,j,k,0);
2050 HypreNodeLap::Int nelems_old = nelems;
2051 ++nelems;
2052
2053 if (gid(i-1,j-1,k) < gidmax) {
2054 cols[nelems] = gid(i-1,j-1,k);
2055 mat[nelems] = sten(i-1,j-1,k,3);
2056 ++nelems;
2057 }
2058
2059 if (gid(i,j-1,k) < gidmax) {
2060 cols[nelems] = gid(i,j-1,k);
2061 mat[nelems] = sten(i,j-1,k,2);
2062 ++nelems;
2063 }
2064
2065 if (gid(i+1,j-1,k) < gidmax) {
2066 cols[nelems] = gid(i+1,j-1,k);
2067 mat[nelems] = sten(i ,j-1,k,3);
2068 ++nelems;
2069 }
2070
2071 if (gid(i-1,j,k) < gidmax) {
2072 cols[nelems] = gid(i-1,j,k);
2073 mat[nelems] = sten(i-1,j,k,1);
2074 ++nelems;
2075 }
2076
2077 if (gid(i+1,j,k) < gidmax) {
2078 cols[nelems] = gid(i+1,j,k);
2079 mat[nelems] = sten(i ,j,k,1);
2080 ++nelems;
2081 }
2082
2083 if (gid(i-1,j+1,k) < gidmax) {
2084 cols[nelems] = gid(i-1,j+1,k);
2085 mat[nelems] = sten(i-1,j ,k,3);
2086 ++nelems;
2087 }
2088
2089 if (gid(i,j+1,k) < gidmax) {
2090 cols[nelems] = gid(i,j+1,k);
2091 mat[nelems] = sten(i,j ,k,2);
2092 ++nelems;
2093 }
2094
2095 if (gid(i+1,j+1,k) < gidmax) {
2096 cols[nelems] = gid(i+1,j+1,k);
2097 mat[nelems] = sten(i ,j ,k,3);
2098 ++nelems;
2099 }
2100
2101 ncols[lid(i,j,k)] = nelems - nelems_old;
2102 }
2103 });
2104}
2105
2106template <typename HypreInt, typename AtomicInt>
2107void mlndlap_fillijmat_aa_cpu (Box const& ndbx,
2108 Array4<AtomicInt const> const& gid,
2109 Array4<int const> const& lid,
2110 HypreInt* ncols, HypreInt* cols,
2111 Real* mat, // NOLINT(readability-non-const-parameter)
2112 Array4<Real const> const& sig,
2113 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2114 Box const& ccdom, bool is_rz) noexcept
2115{
2116 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2117 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2118 Real fxy = facx + facy;
2119 Real f2xmy = Real(2.0)*facx - facy;
2120 Real fmx2y = Real(2.0)*facy - facx;
2121
2122 // Note that ccdom has been grown at periodic boundaries.
2123 const Box& nddom = amrex::surroundingNodes(ccdom);
2124
2125 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2126 HypreInt nelems = 0;
2127 amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2128 {
2129 if (lid(i,j,k) >= 0)
2130 {
2131 Real fp, fm;
2132 if (is_rz) {
2133 fp = facy / static_cast<Real>(2*i+1);
2134 fm = facy / static_cast<Real>(2*i-1);
2135 } else {
2136 fp = fm = Real(0.0);
2137 }
2138
2139 HypreInt nelems_old = nelems;
2140 cols[nelems_old] = gid(i,j,k);
2141 Real m0 = Real(0.);
2142 ++nelems;
2143
2144 if (nddom.contains(i-1,j-1,k)) {
2145 Real tmp = fxy*sig(i-1,j-1,k);
2146 m0 -= tmp;
2147 if ( gid(i-1,j-1,k) < gidmax) {
2148 cols[nelems] = gid(i-1,j-1,k);
2149 mat[nelems] = tmp;
2150 ++nelems;
2151 }
2152 }
2153
2154 if (nddom.contains(i,j-1,k)) {
2155 Real tmp = Real(0.0);
2156 if ( ccdom.contains(i-1,j-1,k)) {
2157 tmp += sig(i-1,j-1,k) * (fmx2y + fm);
2158 }
2159 if ( ccdom.contains(i,j-1,k)) {
2160 tmp += sig(i,j-1,k) * (fmx2y - fp);
2161 }
2162 m0 -= tmp;
2163 if (gid(i,j-1,k) < gidmax) {
2164 cols[nelems] = gid(i,j-1,k);
2165 mat[nelems] = tmp;
2166 ++nelems;
2167 }
2168 }
2169
2170 if (nddom.contains(i+1,j-1,k)) {
2171 Real tmp = fxy*sig(i ,j-1,k);
2172 m0 -= tmp;
2173 if (gid(i+1,j-1,k) < gidmax) {
2174 cols[nelems] = gid(i+1,j-1,k);
2175 mat[nelems] = tmp;
2176 ++nelems;
2177 }
2178 }
2179
2180 if (nddom.contains(i-1,j,k)) {
2181 Real tmp = Real(0.0);
2182 if ( ccdom.contains(i-1,j-1,k)) {
2183 tmp += f2xmy*sig(i-1,j-1,k);
2184 }
2185 if ( ccdom.contains(i-1,j,k)) {
2186 tmp += f2xmy*sig(i-1,j,k);
2187 }
2188 m0 -= tmp;
2189 if (gid(i-1,j,k) < gidmax) {
2190 cols[nelems] = gid(i-1,j,k);
2191 mat[nelems] = tmp;
2192 ++nelems;
2193 }
2194 }
2195
2196 if (nddom.contains(i+1,j,k)) {
2197 Real tmp = Real(0.0);
2198 if ( ccdom.contains(i ,j-1,k)) {
2199 tmp += f2xmy*sig(i ,j-1,k);
2200 }
2201 if ( ccdom.contains(i ,j,k)) {
2202 tmp += f2xmy*sig(i ,j,k);
2203 }
2204 m0 -= tmp;
2205 if (gid(i+1,j,k) < gidmax) {
2206 cols[nelems] = gid(i+1,j,k);
2207 mat[nelems] = tmp;
2208 ++nelems;
2209 }
2210 }
2211
2212 if (nddom.contains(i-1,j+1,k)) {
2213 Real tmp = fxy*sig(i-1,j ,k);
2214 m0 -= tmp;
2215 if (gid(i-1,j+1,k) < gidmax) {
2216 cols[nelems] = gid(i-1,j+1,k);
2217 mat[nelems] = tmp;
2218 ++nelems;
2219 }
2220 }
2221
2222 if (nddom.contains(i,j+1,k)) {
2223 Real tmp = Real(0.0);
2224 if ( ccdom.contains(i-1,j ,k)) {
2225 tmp += sig(i-1,j ,k) * (fmx2y + fm);
2226 }
2227 if ( ccdom.contains(i,j ,k)) {
2228 tmp += sig(i,j ,k) * (fmx2y - fp);
2229 }
2230 m0 -= tmp;
2231 if (gid(i,j+1,k) < gidmax) {
2232 cols[nelems] = gid(i,j+1,k);
2233 mat[nelems] = tmp;
2234 ++nelems;
2235 }
2236 }
2237
2238 if (nddom.contains(i+1,j+1,k)) {
2239 Real tmp = fxy*sig(i ,j ,k);
2240 m0 -= tmp;
2241 if (gid(i+1,j+1,k) < gidmax) {
2242 cols[nelems] = gid(i+1,j+1,k);
2243 mat[nelems] = tmp;
2244 ++nelems;
2245 }
2246 }
2247
2248 mat[nelems_old] = m0;
2249 ncols[lid(i,j,k)] = nelems - nelems_old;
2250 }
2251 });
2252}
2253
2254template <typename HypreInt, typename AtomicInt>
2255void mlndlap_fillijmat_ha_cpu (Box const& ndbx,
2256 Array4<AtomicInt const> const& gid,
2257 Array4<int const> const& lid,
2258 HypreInt* ncols, HypreInt* cols,
2259 Real* mat, // NOLINT(readability-non-const-parameter)
2260 Array4<Real const> const& sx,
2261 Array4<Real const> const& sy,
2262 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2263 Box const& ccdom, bool is_rz) noexcept
2264{
2265 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2266 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2267
2268 // Note that ccdom has been grown at periodic boundaries.
2269 const Box& nddom = amrex::surroundingNodes(ccdom);
2270
2271 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2272 HypreInt nelems = 0;
2273 amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2274 {
2275 if (lid(i,j,k) >= 0)
2276 {
2277 Real fp, fm;
2278 if (is_rz) {
2279 fp = facy / static_cast<Real>(2*i+1);
2280 fm = facy / static_cast<Real>(2*i-1);
2281 } else {
2282 fp = fm = Real(0.0);
2283 }
2284
2285 HypreInt nelems_old = nelems;
2286 cols[nelems_old] = gid(i,j,k);
2287 Real m0 = Real(0.);
2288 ++nelems;
2289
2290 if (nddom.contains(i-1,j-1,k)) {
2291 Real tmp = facx*sx(i-1,j-1,k) + facy*sy(i-1,j-1,k);
2292 m0 -= tmp;
2293 if ( gid(i-1,j-1,k) < gidmax) {
2294 cols[nelems] = gid(i-1,j-1,k);
2295 mat[nelems] = tmp;
2296 ++nelems;
2297 }
2298 }
2299
2300 if (nddom.contains(i,j-1,k)) {
2301 Real tmp = Real(0.0);
2302 if ( ccdom.contains(i-1,j-1,k)) {
2303 tmp += sy(i-1,j-1,k) * (facy * Real(2.0) + fm)
2304 - sx(i-1,j-1,k) * facx;
2305 }
2306 if ( ccdom.contains(i,j-1,k)) {
2307 tmp += sy(i,j-1,k) * (facy * Real(2.0) - fp)
2308 - sx(i,j-1,k) * facx;
2309 }
2310 m0 -= tmp;
2311 if (gid(i,j-1,k) < gidmax) {
2312 cols[nelems] = gid(i,j-1,k);
2313 mat[nelems] = tmp;
2314 ++nelems;
2315 }
2316 }
2317
2318 if (nddom.contains(i+1,j-1,k)) {
2319 Real tmp = facx*sx(i ,j-1,k) + facy*sy(i ,j-1,k);
2320 m0 -= tmp;
2321 if (gid(i+1,j-1,k) < gidmax) {
2322 cols[nelems] = gid(i+1,j-1,k);
2323 mat[nelems] = tmp;
2324 ++nelems;
2325 }
2326 }
2327
2328 if (nddom.contains(i-1,j,k)) {
2329 Real tmp = Real(0.0);
2330 if ( ccdom.contains(i-1,j-1,k)) {
2331 tmp += sx(i-1,j-1,k) * facx*Real(2.0)
2332 - sy(i-1,j-1,k) * facy;
2333 }
2334 if ( ccdom.contains(i-1,j,k)) {
2335 tmp += sx(i-1,j,k) * facx*Real(2.0)
2336 - sy(i-1,j,k) * facy;
2337 }
2338 m0 -= tmp;
2339 if (gid(i-1,j,k) < gidmax) {
2340 cols[nelems] = gid(i-1,j,k);
2341 mat[nelems] = tmp;
2342 ++nelems;
2343 }
2344 }
2345
2346 if (nddom.contains(i+1,j,k)) {
2347 Real tmp = Real(0.0);
2348 if ( ccdom.contains(i ,j-1,k)) {
2349 tmp += sx(i ,j-1,k) * facx*Real(2.0)
2350 - sy(i ,j-1,k) * facy;
2351 }
2352 if ( ccdom.contains(i ,j,k)) {
2353 tmp += sx(i ,j,k) * facx*Real(2.0)
2354 - sy(i ,j,k) * facy;
2355 }
2356 m0 -= tmp;
2357 if (gid(i+1,j,k) < gidmax) {
2358 cols[nelems] = gid(i+1,j,k);
2359 mat[nelems] = tmp;
2360 ++nelems;
2361 }
2362 }
2363
2364 if (nddom.contains(i-1,j+1,k)) {
2365 Real tmp = facx*sx(i-1,j ,k) + facy*sy(i-1,j ,k);
2366 m0 -= tmp;
2367 if (gid(i-1,j+1,k) < gidmax) {
2368 cols[nelems] = gid(i-1,j+1,k);
2369 mat[nelems] = tmp;
2370 ++nelems;
2371 }
2372 }
2373
2374 if (nddom.contains(i,j+1,k)) {
2375 Real tmp = Real(0.0);
2376 if ( ccdom.contains(i-1,j ,k)) {
2377 tmp += sy(i-1,j ,k) * (facy*Real(2.0) + fm)
2378 - sx(i-1,j ,k) * facx;
2379 }
2380 if ( ccdom.contains(i,j ,k)) {
2381 tmp += sy(i,j ,k) * (facy*Real(2.0) - fp)
2382 - sx(i,j ,k) * facx;
2383 }
2384 m0 -= tmp;
2385 if (gid(i,j+1,k) < gidmax) {
2386 cols[nelems] = gid(i,j+1,k);
2387 mat[nelems] = tmp;
2388 ++nelems;
2389 }
2390 }
2391
2392 if (nddom.contains(i+1,j+1,k)) {
2393 Real tmp = facx*sx(i ,j ,k) + facy*sy(i ,j ,k);
2394 m0 -= tmp;
2395 if (gid(i+1,j+1,k) < gidmax) {
2396 cols[nelems] = gid(i+1,j+1,k);
2397 mat[nelems] = tmp;
2398 ++nelems;
2399 }
2400 }
2401
2402 mat[nelems_old] = m0;
2403 ncols[lid(i,j,k)] = nelems - nelems_old;
2404 }
2405 });
2406}
2407
2408template <typename HypreInt, typename AtomicInt>
2409void mlndlap_fillijmat_cs_cpu (Box const& ndbx,
2410 Array4<AtomicInt const> const& gid,
2411 Array4<int const> const& lid,
2412 HypreInt* ncols, HypreInt* cols,
2413 Real* mat, // NOLINT(readability-non-const-parameter)
2414 Real sig,
2415 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2416 Box const& ccdom, bool is_rz) noexcept
2417{
2418 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0] * sig;
2419 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1] * sig;
2420 Real fxy = facx + facy;
2421 Real f2xmy = Real(2.0)*facx - facy;
2422 Real fmx2y = Real(2.0)*facy - facx;
2423
2424 // Note that ccdom has been grown at periodic boundaries.
2425 const Box& nddom = amrex::surroundingNodes(ccdom);
2426
2427 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2428 HypreInt nelems = 0;
2429 amrex::LoopOnCpu(ndbx, [&] (int i, int j, int k) noexcept
2430 {
2431 if (lid(i,j,k) >= 0)
2432 {
2433 Real fp, fm;
2434 if (is_rz) {
2435 fp = facy / static_cast<Real>(2*i+1);
2436 fm = facy / static_cast<Real>(2*i-1);
2437 } else {
2438 fp = fm = Real(0.0);
2439 }
2440
2441 HypreInt nelems_old = nelems;
2442 cols[nelems_old] = gid(i,j,k);
2443 Real m0 = Real(0.);
2444 ++nelems;
2445
2446 if (nddom.contains(i-1,j-1,k)) {
2447 Real tmp = fxy;
2448 m0 -= tmp;
2449 if ( gid(i-1,j-1,k) < gidmax) {
2450 cols[nelems] = gid(i-1,j-1,k);
2451 mat[nelems] = tmp;
2452 ++nelems;
2453 }
2454 }
2455
2456 if (nddom.contains(i,j-1,k)) {
2457 Real tmp = Real(0.0);
2458 if ( ccdom.contains(i-1,j-1,k)) {
2459 tmp += fmx2y + fm;
2460 }
2461 if ( ccdom.contains(i,j-1,k)) {
2462 tmp += fmx2y - fp;
2463 }
2464 m0 -= tmp;
2465 if (gid(i,j-1,k) < gidmax) {
2466 cols[nelems] = gid(i,j-1,k);
2467 mat[nelems] = tmp;
2468 ++nelems;
2469 }
2470 }
2471
2472 if (nddom.contains(i+1,j-1,k)) {
2473 Real tmp = fxy;
2474 m0 -= tmp;
2475 if (gid(i+1,j-1,k) < gidmax) {
2476 cols[nelems] = gid(i+1,j-1,k);
2477 mat[nelems] = tmp;
2478 ++nelems;
2479 }
2480 }
2481
2482 if (nddom.contains(i-1,j,k)) {
2483 Real tmp = Real(0.0);
2484 if ( ccdom.contains(i-1,j-1,k)) {
2485 tmp += f2xmy;
2486 }
2487 if ( ccdom.contains(i-1,j,k)) {
2488 tmp += f2xmy;
2489 }
2490 m0 -= tmp;
2491 if (gid(i-1,j,k) < gidmax) {
2492 cols[nelems] = gid(i-1,j,k);
2493 mat[nelems] = tmp;
2494 ++nelems;
2495 }
2496 }
2497
2498 if (nddom.contains(i+1,j,k)) {
2499 Real tmp = Real(0.0);
2500 if ( ccdom.contains(i ,j-1,k)) {
2501 tmp += f2xmy;
2502 }
2503 if ( ccdom.contains(i ,j,k)) {
2504 tmp += f2xmy;
2505 }
2506 m0 -= tmp;
2507 if (gid(i+1,j,k) < gidmax) {
2508 cols[nelems] = gid(i+1,j,k);
2509 mat[nelems] = tmp;
2510 ++nelems;
2511 }
2512 }
2513
2514 if (nddom.contains(i-1,j+1,k)) {
2515 Real tmp = fxy;
2516 m0 -= tmp;
2517 if (gid(i-1,j+1,k) < gidmax) {
2518 cols[nelems] = gid(i-1,j+1,k);
2519 mat[nelems] = tmp;
2520 ++nelems;
2521 }
2522 }
2523
2524 if (nddom.contains(i,j+1,k)) {
2525 Real tmp = Real(0.0);
2526 if ( ccdom.contains(i-1,j ,k)) {
2527 tmp += fmx2y + fm;
2528 }
2529 if ( ccdom.contains(i,j ,k)) {
2530 tmp += fmx2y - fp;
2531 }
2532 m0 -= tmp;
2533 if (gid(i,j+1,k) < gidmax) {
2534 cols[nelems] = gid(i,j+1,k);
2535 mat[nelems] = tmp;
2536 ++nelems;
2537 }
2538 }
2539
2540 if (nddom.contains(i+1,j+1,k)) {
2541 Real tmp = fxy;
2542 m0 -= tmp;
2543 if (gid(i+1,j+1,k) < gidmax) {
2544 cols[nelems] = gid(i+1,j+1,k);
2545 mat[nelems] = tmp;
2546 ++nelems;
2547 }
2548 }
2549
2550 mat[nelems_old] = m0;
2551 ncols[lid(i,j,k)] = nelems - nelems_old;
2552 }
2553 });
2554}
2555
2556#ifdef AMREX_USE_GPU
2557
2558template <typename HypreInt, typename AtomicInt>
2559AMREX_GPU_DEVICE AMREX_FORCE_INLINE
2560void mlndlap_fillijmat_sten_gpu (const int ps, const int i, const int j, const int k,
2561 const int offset,
2562 Array4<AtomicInt const> const& gid,
2563 Array4<int const> const& lid,
2564 HypreInt* ncols, HypreInt* cols,
2565 Real* mat, // NOLINT(readability-non-const-parameter)
2566 Array4<Real const> const& sten) noexcept
2567{
2568 if (lid(i,j,k) >= 0)
2569 {
2570 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2571 int nelems = 0;
2572
2573 if (offset == 1 || offset == 0) {
2574 if (gid(i-1,j-1,k) < gidmax) {
2575 if (offset != 0) {
2576 cols[ps] = gid(i-1,j-1,k);
2577 mat[ps] = sten(i-1,j-1,k,3);
2578 }
2579 ++nelems;
2580 }
2581 if (offset != 0) { return; }
2582 }
2583
2584 if (offset == 2 || offset == 0) {
2585 if (gid(i,j-1,k) < gidmax) {
2586 if (offset != 0) {
2587 cols[ps] = gid(i,j-1,k);
2588 mat[ps] = sten(i,j-1,k,2);
2589 }
2590 ++nelems;
2591 }
2592 if (offset != 0) { return; }
2593 }
2594
2595 if (offset == 3 || offset == 0) {
2596 if (gid(i+1,j-1,k) < gidmax) {
2597 if (offset != 0) {
2598 cols[ps] = gid(i+1,j-1,k);
2599 mat[ps] = sten(i ,j-1,k,3);
2600 }
2601 ++nelems;
2602 }
2603 if (offset != 0) { return; }
2604 }
2605
2606 if (offset == 4 || offset == 0) {
2607 if (gid(i-1,j,k) < gidmax) {
2608 if (offset != 0) {
2609 cols[ps] = gid(i-1,j,k);
2610 mat[ps] = sten(i-1,j,k,1);
2611 }
2612 ++nelems;
2613 }
2614 if (offset != 0) { return; }
2615 }
2616
2617 if (offset == 5 || offset == 0) {
2618 if (gid(i+1,j,k) < gidmax) {
2619 if (offset != 0) {
2620 cols[ps] = gid(i+1,j,k);
2621 mat[ps] = sten(i ,j,k,1);
2622 }
2623 ++nelems;
2624 }
2625 if (offset != 0) { return; }
2626 }
2627
2628 if (offset == 6 || offset == 0) {
2629 if (gid(i-1,j+1,k) < gidmax) {
2630 if (offset != 0) {
2631 cols[ps] = gid(i-1,j+1,k);
2632 mat[ps] = sten(i-1,j ,k,3);
2633 }
2634 ++nelems;
2635 }
2636 if (offset != 0) { return; }
2637 }
2638
2639 if (offset == 7 || offset == 0) {
2640 if (gid(i,j+1,k) < gidmax) {
2641 if (offset != 0) {
2642 cols[ps] = gid(i,j+1,k);
2643 mat[ps] = sten(i,j ,k,2);
2644 }
2645 ++nelems;
2646 }
2647 if (offset != 0) { return; }
2648 }
2649
2650 if (offset == 8 || offset == 0) {
2651 if (gid(i+1,j+1,k) < gidmax) {
2652 if (offset != 0) {
2653 cols[ps] = gid(i+1,j+1,k);
2654 mat[ps] = sten(i ,j ,k,3);
2655 }
2656 ++nelems;
2657 }
2658 if (offset != 0) { return; }
2659 }
2660
2661 // Only offset == 0 could get this far.
2662 cols[ps] = gid(i,j,k);
2663 mat[ps] = sten(i,j,k,0);
2664 ncols[lid(i,j,k)] = nelems+1;
2665 }
2666}
2667
2668template <typename HypreInt, typename AtomicInt>
2669AMREX_GPU_DEVICE AMREX_FORCE_INLINE
2670void mlndlap_fillijmat_aa_gpu (const int ps, const int i, const int j, const int k,
2671 const int offset,
2672 Box const& ndbx, Array4<AtomicInt const> const& gid,
2673 Array4<int const> const& lid,
2674 HypreInt* ncols, HypreInt* cols,
2675 Real* mat, // NOLINT(readability-non-const-parameter)
2676 Array4<Real const> const& sig,
2677 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2678 Box const& ccdom, bool is_rz) noexcept
2679{
2680 if (lid(i,j,k) >= 0)
2681 {
2682 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2683 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2684 Real fxy = facx + facy;
2685 Real f2xmy = Real(2.0)*facx - facy;
2686 Real fmx2y = Real(2.0)*facy - facx;
2687
2688 Real fp, fm;
2689 if (is_rz) {
2690 fp = facy / static_cast<Real>(2*i+1);
2691 fm = facy / static_cast<Real>(2*i-1);
2692 } else {
2693 fp = fm = Real(0.0);
2694 }
2695
2696 const Box& nddom = amrex::surroundingNodes(ccdom);
2697
2698 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2699 int nelems = 0;
2700 Real m0 = Real(0.);
2701
2702 if (offset == 1 || offset == 0) {
2703 if (nddom.contains(i-1,j-1,k)) {
2704 Real tmp = fxy*sig(i-1,j-1,k);
2705 m0 -= tmp;
2706 if (gid(i-1,j-1,k) < gidmax) {
2707 if (offset != 0) {
2708 cols[ps] = gid(i-1,j-1,k);
2709 mat[ps] = tmp;
2710 }
2711 ++nelems;
2712 }
2713 }
2714 if (offset != 0) { return; }
2715 }
2716
2717 if (offset == 2 || offset == 0) {
2718 if (nddom.contains(i,j-1,k)) {
2719 Real tmp = Real(0.0);
2720 if ( ccdom.contains(i-1,j-1,k)) {
2721 tmp += sig(i-1,j-1,k) * (fmx2y + fm);
2722 }
2723 if ( ccdom.contains(i,j-1,k)) {
2724 tmp += sig(i,j-1,k) * (fmx2y - fp);
2725 }
2726 m0 -= tmp;
2727 if (gid(i,j-1,k) < gidmax) {
2728 if (offset != 0) {
2729 cols[ps] = gid(i,j-1,k);
2730 mat[ps] = tmp;
2731 }
2732 ++nelems;
2733 }
2734 }
2735 if (offset != 0) { return; }
2736 }
2737
2738 if (offset == 3 || offset == 0) {
2739 if (nddom.contains(i+1,j-1,k)) {
2740 Real tmp = fxy*sig(i ,j-1,k);
2741 m0 -= tmp;
2742 if (gid(i+1,j-1,k) < gidmax) {
2743 if (offset != 0) {
2744 cols[ps] = gid(i+1,j-1,k);
2745 mat[ps] = tmp;
2746 }
2747 ++nelems;
2748 }
2749 }
2750 if (offset != 0) { return; }
2751 }
2752
2753 if (offset == 4 || offset == 0) {
2754 if (nddom.contains(i-1,j,k)) {
2755 Real tmp = Real(0.0);
2756 if ( ccdom.contains(i-1,j-1,k)) {
2757 tmp += f2xmy*sig(i-1,j-1,k);
2758 }
2759 if ( ccdom.contains(i-1,j,k)) {
2760 tmp += f2xmy*sig(i-1,j,k);
2761 }
2762 m0 -= tmp;
2763 if (gid(i-1,j,k) < gidmax) {
2764 if (offset != 0) {
2765 cols[ps] = gid(i-1,j,k);
2766 mat[ps] = tmp;
2767 }
2768 ++nelems;
2769 }
2770 }
2771 if (offset != 0) { return; }
2772 }
2773
2774 if (offset == 5 || offset == 0) {
2775 if (nddom.contains(i+1,j,k)) {
2776 Real tmp = Real(0.0);
2777 if ( ccdom.contains(i ,j-1,k)) {
2778 tmp += f2xmy*sig(i ,j-1,k);
2779 }
2780 if ( ccdom.contains(i ,j,k)) {
2781 tmp += f2xmy*sig(i ,j,k);
2782 }
2783 m0 -= tmp;
2784 if (gid(i+1,j,k) < gidmax) {
2785 if (offset != 0) {
2786 cols[ps] = gid(i+1,j,k);
2787 mat[ps] = tmp;
2788 }
2789 ++nelems;
2790 }
2791 }
2792 if (offset != 0) { return; }
2793 }
2794
2795 if (offset == 6 || offset == 0) {
2796 if (nddom.contains(i-1,j+1,k)) {
2797 Real tmp = fxy*sig(i-1,j ,k);
2798 m0 -= tmp;
2799 if (gid(i-1,j+1,k) < gidmax) {
2800 if (offset != 0) {
2801 cols[ps] = gid(i-1,j+1,k);
2802 mat[ps] = tmp;
2803 }
2804 ++nelems;
2805 }
2806 }
2807 if (offset != 0) { return; }
2808 }
2809
2810 if (offset == 7 || offset == 0) {
2811 if (nddom.contains(i,j+1,k)) {
2812 Real tmp = Real(0.0);
2813 if ( ccdom.contains(i-1,j ,k)) {
2814 tmp += sig(i-1,j ,k) * (fmx2y + fm);
2815 }
2816 if ( ccdom.contains(i,j ,k)) {
2817 tmp += sig(i,j ,k) * (fmx2y - fp);
2818 }
2819 m0 -= tmp;
2820 if (gid(i,j+1,k) < gidmax) {
2821 if (offset != 0) {
2822 cols[ps] = gid(i,j+1,k);
2823 mat[ps] = tmp;
2824 }
2825 ++nelems;
2826 }
2827 }
2828 if (offset != 0) { return; }
2829 }
2830
2831 if (offset == 8 || offset == 0) {
2832 if (nddom.contains(i+1,j+1,k)) {
2833 Real tmp = fxy*sig(i ,j ,k);
2834 m0 -= tmp;
2835 if (gid(i+1,j+1,k) < gidmax) {
2836 if (offset != 0) {
2837 cols[ps] = gid(i+1,j+1,k);
2838 mat[ps] = tmp;
2839 }
2840 ++nelems;
2841 }
2842 }
2843 if (offset != 0) { return; }
2844 }
2845
2846 // Only offset == 0 could get this far.
2847 cols[ps] = gid(i,j,k);
2848 mat[ps] = m0;
2849 ncols[lid(i,j,k)] = nelems+1;
2850 }
2851}
2852
2853template <typename HypreInt, typename AtomicInt>
2854AMREX_GPU_DEVICE AMREX_FORCE_INLINE
2855void mlndlap_fillijmat_ha_gpu (const int ps, const int i, const int j, const int k,
2856 const int offset,
2857 Box const& ndbx, Array4<AtomicInt const> const& gid,
2858 Array4<int const> const& lid,
2859 HypreInt* ncols, HypreInt* cols,
2860 Real* mat, // NOLINT(readability-non-const-parameter)
2861 Array4<Real const> const& sx,
2862 Array4<Real const> const& sy,
2863 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
2864 Box const& ccdom, bool is_rz) noexcept
2865{
2866 if (lid(i,j,k) >= 0)
2867 {
2868 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
2869 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
2870
2871 Real fp, fm;
2872 if (is_rz) {
2873 fp = facy / static_cast<Real>(2*i+1);
2874 fm = facy / static_cast<Real>(2*i-1);
2875 } else {
2876 fp = fm = Real(0.0);
2877 }
2878
2879 const Box& nddom = amrex::surroundingNodes(ccdom);
2880
2881 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
2882 int nelems = 0;
2883 Real m0 = Real(0.);
2884
2885 if (offset == 1 || offset == 0) {
2886 if (nddom.contains(i-1,j-1,k)) {
2887 Real tmp = facx*sx(i-1,j-1,k) + facy*sy(i-1,j-1,k);
2888 m0 -= tmp;
2889 if ( gid(i-1,j-1,k) < gidmax) {
2890 if (offset != 0) {
2891 cols[ps] = gid(i-1,j-1,k);
2892 mat[ps] = tmp;
2893 }
2894 ++nelems;
2895 }
2896 }
2897 if (offset != 0) { return; }
2898 }
2899
2900 if (offset == 2 || offset == 0) {
2901 if (nddom.contains(i,j-1,k)) {
2902 Real tmp = Real(0.0);
2903 if ( ccdom.contains(i-1,j-1,k)) {
2904 tmp += sy(i-1,j-1,k) * (facy * Real(2.0) + fm)
2905 - sx(i-1,j-1,k) * facx;
2906 }
2907 if ( ccdom.contains(i,j-1,k)) {
2908 tmp += sy(i,j-1,k) * (facy * Real(2.0) - fp)
2909 - sx(i,j-1,k) * facx;
2910 }
2911 m0 -= tmp;
2912 if (gid(i,j-1,k) < gidmax) {
2913 if (offset != 0) {
2914 cols[ps] = gid(i,j-1,k);
2915 mat[ps] = tmp;
2916 }
2917 ++nelems;
2918 }
2919 }
2920 if (offset != 0) { return; }
2921 }
2922
2923 if (offset == 3 || offset == 0) {
2924 if (nddom.contains(i+1,j-1,k)) {
2925 Real tmp = facx*sx(i ,j-1,k) + facy*sy(i ,j-1,k);
2926 m0 -= tmp;
2927 if (gid(i+1,j-1,k) < gidmax) {
2928 if (offset != 0) {
2929 cols[ps] = gid(i+1,j-1,k);
2930 mat[ps] = tmp;
2931 }
2932 ++nelems;
2933 }
2934 }
2935 if (offset != 0) { return; }
2936 }
2937
2938 if (offset == 4 || offset == 0) {
2939 if (nddom.contains(i-1,j,k)) {
2940 Real tmp = Real(0.0);
2941 if ( ccdom.contains(i-1,j-1,k)) {
2942 tmp += sx(i-1,j-1,k) * facx*Real(2.0)
2943 - sy(i-1,j-1,k) * facy;
2944 }
2945 if ( ccdom.contains(i-1,j,k)) {
2946 tmp += sx(i-1,j,k) * facx*Real(2.0)
2947 - sy(i-1,j,k) * facy;
2948 }
2949 m0 -= tmp;
2950 if (gid(i-1,j,k) < gidmax) {
2951 if (offset != 0) {
2952 cols[ps] = gid(i-1,j,k);
2953 mat[ps] = tmp;
2954 }
2955 ++nelems;
2956 }
2957 }
2958 if (offset != 0) { return; }
2959 }
2960
2961 if (offset == 5 || offset == 0) {
2962 if (nddom.contains(i+1,j,k)) {
2963 Real tmp = Real(0.0);
2964 if ( ccdom.contains(i ,j-1,k)) {
2965 tmp += sx(i ,j-1,k) * facx*Real(2.0)
2966 - sy(i ,j-1,k) * facy;
2967 }
2968 if ( ccdom.contains(i ,j,k)) {
2969 tmp += sx(i ,j,k) * facx*Real(2.0)
2970 - sy(i ,j,k) * facy;
2971 }
2972 m0 -= tmp;
2973 if (gid(i+1,j,k) < gidmax) {
2974 if (offset != 0) {
2975 cols[ps] = gid(i+1,j,k);
2976 mat[ps] = tmp;
2977 }
2978 ++nelems;
2979 }
2980 }
2981 if (offset != 0) { return; }
2982 }
2983
2984 if (offset == 6 || offset == 0) {
2985 if (nddom.contains(i-1,j+1,k)) {
2986 Real tmp = facx*sx(i-1,j ,k) + facy*sy(i-1,j ,k);
2987 m0 -= tmp;
2988 if (gid(i-1,j+1,k) < gidmax) {
2989 if (offset != 0) {
2990 cols[ps] = gid(i-1,j+1,k);
2991 mat[ps] = tmp;
2992 }
2993 ++nelems;
2994 }
2995 }
2996 if (offset != 0) { return; }
2997 }
2998
2999 if (offset == 7 || offset == 0) {
3000 if (nddom.contains(i,j+1,k)) {
3001 Real tmp = Real(0.0);
3002 if ( ccdom.contains(i-1,j ,k)) {
3003 tmp += sy(i-1,j ,k) * (facy*Real(2.0) + fm)
3004 - sx(i-1,j ,k) * facx;
3005 }
3006 if ( ccdom.contains(i,j ,k)) {
3007 tmp += sy(i,j ,k) * (facy*Real(2.0) - fp)
3008 - sx(i,j ,k) * facx;
3009 }
3010 m0 -= tmp;
3011 if (gid(i,j+1,k) < gidmax) {
3012 if (offset != 0) {
3013 cols[ps] = gid(i,j+1,k);
3014 mat[ps] = tmp;
3015 }
3016 ++nelems;
3017 }
3018 }
3019 if (offset != 0) { return; }
3020 }
3021
3022 if (offset == 8 || offset == 0) {
3023 if (nddom.contains(i+1,j+1,k)) {
3024 Real tmp = facx*sx(i ,j ,k) + facy*sy(i ,j ,k);
3025 m0 -= tmp;
3026 if (gid(i+1,j+1,k) < gidmax) {
3027 if (offset != 0) {
3028 cols[ps] = gid(i+1,j+1,k);
3029 mat[ps] = tmp;
3030 }
3031 ++nelems;
3032 }
3033 }
3034 if (offset != 0) { return; }
3035 }
3036
3037 // Only offset == 0 could get this far.
3038 cols[ps] = gid(i,j,k);
3039 mat[ps] = m0;
3040 ncols[lid(i,j,k)] = nelems+1;
3041 }
3042}
3043
3044template <typename HypreInt, typename AtomicInt>
3045AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3046void mlndlap_fillijmat_cs_gpu (const int ps, const int i, const int j, const int k,
3047 const int offset,
3048 Box const& ndbx, Array4<AtomicInt const> const& gid,
3049 Array4<int const> const& lid,
3050 HypreInt* ncols, HypreInt* cols,
3051 Real* mat, // NOLINT(readability-non-const-parameter)
3052 Real sig, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
3053 Box const& ccdom, bool is_rz) noexcept
3054{
3055 if (lid(i,j,k) >= 0)
3056 {
3057 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0] * sig;
3058 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1] * sig;
3059 Real fxy = facx + facy;
3060 Real f2xmy = Real(2.0)*facx - facy;
3061 Real fmx2y = Real(2.0)*facy - facx;
3062
3063 Real fp, fm;
3064 if (is_rz) {
3065 fp = facy / static_cast<Real>(2*i+1);
3066 fm = facy / static_cast<Real>(2*i-1);
3067 } else {
3068 fp = fm = Real(0.0);
3069 }
3070
3071 // Note that nddom has been grown at periodic boundaries.
3072 const Box& nddom = amrex::surroundingNodes(ccdom);
3073
3074 constexpr auto gidmax = std::numeric_limits<AtomicInt>::max();
3075 int nelems = 0;
3076 Real m0 = Real(0.);
3077
3078 if (offset == 1 || offset == 0) {
3079 if (nddom.contains(i-1,j-1,k)) {
3080 Real tmp = fxy;
3081 m0 -= tmp;
3082 if ( gid(i-1,j-1,k) < gidmax) {
3083 if (offset != 0) {
3084 cols[ps] = gid(i-1,j-1,k);
3085 mat[ps] = tmp;
3086 }
3087 ++nelems;
3088 }
3089 }
3090 if (offset != 0) { return; }
3091 }
3092
3093 if (offset == 2 || offset == 0) {
3094 if (nddom.contains(i,j-1,k)) {
3095 Real tmp = Real(0.0);
3096 if ( ccdom.contains(i-1,j-1,k)) {
3097 tmp += fmx2y + fm;
3098 }
3099 if ( ccdom.contains(i,j-1,k)) {
3100 tmp += fmx2y - fp;
3101 }
3102 m0 -= tmp;
3103 if (gid(i,j-1,k) < gidmax) {
3104 if (offset != 0) {
3105 cols[ps] = gid(i,j-1,k);
3106 mat[ps] = tmp;
3107 }
3108 ++nelems;
3109 }
3110 }
3111 if (offset != 0) { return; }
3112 }
3113
3114 if (offset == 3 || offset == 0) {
3115 if (nddom.contains(i+1,j-1,k)) {
3116 Real tmp = fxy;
3117 m0 -= tmp;
3118 if (gid(i+1,j-1,k) < gidmax) {
3119 if (offset != 0) {
3120 cols[ps] = gid(i+1,j-1,k);
3121 mat[ps] = tmp;
3122 }
3123 ++nelems;
3124 }
3125 }
3126 if (offset != 0) { return; }
3127 }
3128
3129 if (offset == 4 || offset == 0) {
3130 if (nddom.contains(i-1,j,k)) {
3131 Real tmp = Real(0.0);
3132 if ( ccdom.contains(i-1,j-1,k)) {
3133 tmp += f2xmy;
3134 }
3135 if ( ccdom.contains(i-1,j,k)) {
3136 tmp += f2xmy;
3137 }
3138 m0 -= tmp;
3139 if (gid(i-1,j,k) < gidmax) {
3140 if (offset != 0) {
3141 cols[ps] = gid(i-1,j,k);
3142 mat[ps] = tmp;
3143 }
3144 ++nelems;
3145 }
3146 }
3147 if (offset != 0) { return; }
3148 }
3149
3150 if (offset == 5 || offset == 0) {
3151 if (nddom.contains(i+1,j,k)) {
3152 Real tmp = Real(0.0);
3153 if ( ccdom.contains(i ,j-1,k)) {
3154 tmp += f2xmy;
3155 }
3156 if ( ccdom.contains(i ,j,k)) {
3157 tmp += f2xmy;
3158 }
3159 m0 -= tmp;
3160 if (gid(i+1,j,k) < gidmax) {
3161 if (offset != 0) {
3162 cols[ps] = gid(i+1,j,k);
3163 mat[ps] = tmp;
3164 }
3165 ++nelems;
3166 }
3167 }
3168 if (offset != 0) { return; }
3169 }
3170
3171 if (offset == 6 || offset == 0) {
3172 if (nddom.contains(i-1,j+1,k)) {
3173 Real tmp = fxy;
3174 m0 -= tmp;
3175 if (gid(i-1,j+1,k) < gidmax) {
3176 if (offset != 0) {
3177 cols[ps] = gid(i-1,j+1,k);
3178 mat[ps] = tmp;
3179 }
3180 ++nelems;
3181 }
3182 }
3183 if (offset != 0) { return; }
3184 }
3185
3186 if (offset == 7 || offset == 0) {
3187 if (nddom.contains(i,j+1,k)) {
3188 Real tmp = Real(0.0);
3189 if ( ccdom.contains(i-1,j ,k)) {
3190 tmp += fmx2y + fm;
3191 }
3192 if ( ccdom.contains(i,j ,k)) {
3193 tmp += fmx2y - fp;
3194 }
3195 m0 -= tmp;
3196 if (gid(i,j+1,k) < gidmax) {
3197 if (offset != 0) {
3198 cols[ps] = gid(i,j+1,k);
3199 mat[ps] = tmp;
3200 }
3201 ++nelems;
3202 }
3203 }
3204 if (offset != 0) { return; }
3205 }
3206
3207 if (offset == 8 || offset == 0) {
3208 if (nddom.contains(i+1,j+1,k)) {
3209 Real tmp = fxy;
3210 m0 -= tmp;
3211 if (gid(i+1,j+1,k) < gidmax) {
3212 if (offset != 0) {
3213 cols[ps] = gid(i+1,j+1,k);
3214 mat[ps] = tmp;
3215 }
3216 ++nelems;
3217 }
3218 }
3219 if (offset != 0) { return; }
3220 }
3221
3222 // Only offset == 0 could get this far.
3223 cols[ps] = gid(i,j,k);
3224 mat[ps] = m0;
3225 ncols[lid(i,j,k)] = nelems+1;
3226 }
3227}
3228
3229#endif
3230
3231#endif
3232
3233AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3234int mlndlap_color (int i, int j, int)
3235{
3236 return (i%2) + (j%2)*2;
3237}
3238
3239AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3240void mlndlap_gscolor_ha (int i, int j, int k, Array4<Real> const& sol,
3241 Array4<Real const> const& rhs, Array4<Real const> const& sx,
3242 Array4<Real const> const& sy, Array4<int const> const& msk,
3243 GpuArray<Real,AMREX_SPACEDIM> const& dxinv, int color,
3244 bool is_rz) noexcept
3245{
3246 if (mlndlap_color(i,j,k) == color) {
3247 if (msk(i,j,k)) {
3248 sol(i,j,k) = Real(0.0);
3249 } else {
3250 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3251 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3252
3253 Real s0 = Real(-2.0)*(facx*(sx(i-1,j-1,k)+sx(i,j-1,k)+sx(i-1,j,k)+sx(i,j,k))
3254 +facy*(sy(i-1,j-1,k)+sy(i,j-1,k)+sy(i-1,j,k)+sy(i,j,k)));
3255
3256 Real Ax = sol(i-1,j-1,k)*(facx*sx(i-1,j-1,k)+facy*sy(i-1,j-1,k))
3257 + sol(i+1,j-1,k)*(facx*sx(i ,j-1,k)+facy*sy(i ,j-1,k))
3258 + sol(i-1,j+1,k)*(facx*sx(i-1,j ,k)+facy*sy(i-1,j ,k))
3259 + sol(i+1,j+1,k)*(facx*sx(i ,j ,k)+facy*sy(i ,j ,k))
3260 + sol(i-1,j,k)*(Real(2.0)*facx*(sx(i-1,j-1,k)+sx(i-1,j,k))
3261 - facy*(sy(i-1,j-1,k)+sy(i-1,j,k)))
3262 + sol(i+1,j,k)*(Real(2.0)*facx*(sx(i ,j-1,k)+sx(i ,j,k))
3263 - facy*(sy(i ,j-1,k)+sy(i ,j,k)))
3264 + sol(i,j-1,k)*( -facx*(sx(i-1,j-1,k)+sx(i,j-1,k))
3265 +Real(2.0)*facy*(sy(i-1,j-1,k)+sy(i,j-1,k)))
3266 + sol(i,j+1,k)*( -facx*(sx(i-1,j ,k)+sx(i,j ,k))
3267 +Real(2.0)*facy*(sy(i-1,j ,k)+sy(i,j ,k)))
3268 + sol(i,j,k)*s0;
3269
3270 if (is_rz) {
3271 Real fp = facy / static_cast<Real>(2*i+1);
3272 Real fm = facy / static_cast<Real>(2*i-1);
3273 Real frzlo = fm*sy(i-1,j-1,k)-fp*sy(i,j-1,k);
3274 Real frzhi = fm*sy(i-1,j ,k)-fp*sy(i,j ,k);
3275 s0 += - frzhi - frzlo;
3276 Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
3277 + frzlo*(sol(i,j-1,k)-sol(i,j,k));
3278 }
3279
3280 sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
3281 }
3282 }
3283}
3284
3285AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3286void mlndlap_gscolor_aa (int i, int j, int k, Array4<Real> const& sol,
3287 Array4<Real const> const& rhs, Array4<Real const> const& sig,
3288 Array4<int const> const& msk,
3289 GpuArray<Real,AMREX_SPACEDIM> const& dxinv, int color,
3290 bool is_rz) noexcept
3291{
3292 if (mlndlap_color(i,j,k) == color) {
3293 if (msk(i,j,k)) {
3294 sol(i,j,k) = Real(0.0);
3295 } else {
3296 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3297 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3298 Real fxy = facx + facy;
3299 Real f2xmy = Real(2.0)*facx - facy;
3300 Real fmx2y = Real(2.0)*facy - facx;
3301
3302 Real s0 = (-Real(2.0))*fxy*(sig(i-1,j-1,k)+sig(i,j-1,k)+sig(i-1,j,k)+sig(i,j,k));
3303 Real Ax = sol(i-1,j-1,k)*fxy*sig(i-1,j-1,k)
3304 + sol(i+1,j-1,k)*fxy*sig(i ,j-1,k)
3305 + sol(i-1,j+1,k)*fxy*sig(i-1,j ,k)
3306 + sol(i+1,j+1,k)*fxy*sig(i ,j ,k)
3307 + sol(i-1,j,k)*f2xmy*(sig(i-1,j-1,k)+sig(i-1,j,k))
3308 + sol(i+1,j,k)*f2xmy*(sig(i ,j-1,k)+sig(i ,j,k))
3309 + sol(i,j-1,k)*fmx2y*(sig(i-1,j-1,k)+sig(i,j-1,k))
3310 + sol(i,j+1,k)*fmx2y*(sig(i-1,j ,k)+sig(i,j ,k))
3311 + sol(i,j,k)*s0;
3312
3313 if (is_rz) {
3314 Real fp = facy / static_cast<Real>(2*i+1);
3315 Real fm = facy / static_cast<Real>(2*i-1);
3316 Real frzlo = fm*sig(i-1,j-1,k)-fp*sig(i,j-1,k);
3317 Real frzhi = fm*sig(i-1,j ,k)-fp*sig(i,j ,k);
3318 s0 += - frzhi - frzlo;
3319 Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
3320 + frzlo*(sol(i,j-1,k)-sol(i,j,k));
3321 }
3322
3323 sol(i,j,k) += (rhs(i,j,k) - Ax) / s0;
3324 }
3325 }
3326}
3327
3328AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3329void mlndlap_gscolor_c (int i, int j, int k, Array4<Real> const& sol,
3330 Array4<Real const> const& rhs, Real sig,
3331 Array4<int const> const& msk,
3332 GpuArray<Real,AMREX_SPACEDIM> const& dxinv, int color,
3333 bool is_rz) noexcept
3334{
3335 if (mlndlap_color(i,j,k) == color) {
3336 if (msk(i,j,k)) {
3337 sol(i,j,k) = Real(0.0);
3338 } else {
3339 Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
3340 Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
3341 Real fxy = facx + facy;
3342 Real f2xmy = Real(2.0)*facx - facy;
3343 Real fmx2y = Real(2.0)*facy - facx;
3344
3345 Real s0 = (-Real(2.0))*fxy*Real(4.);
3346 Real Ax = sol(i-1,j-1,k)*fxy
3347 + sol(i+1,j-1,k)*fxy
3348 + sol(i-1,j+1,k)*fxy
3349 + sol(i+1,j+1,k)*fxy
3350 + sol(i-1,j,k)*f2xmy*Real(2.)
3351 + sol(i+1,j,k)*f2xmy*Real(2.)
3352 + sol(i,j-1,k)*fmx2y*Real(2.)
3353 + sol(i,j+1,k)*fmx2y*Real(2.)
3354 + sol(i,j,k)*s0;
3355
3356 if (is_rz) {
3357 Real fp = facy / static_cast<Real>(2*i+1);
3358 Real fm = facy / static_cast<Real>(2*i-1);
3359 Real frzlo = fm-fp;
3360 Real frzhi = fm-fp;
3361 s0 += - frzhi - frzlo;
3362 Ax += frzhi*(sol(i,j+1,k)-sol(i,j,k))
3363 + frzlo*(sol(i,j-1,k)-sol(i,j,k));
3364 }
3365
3366 sol(i,j,k) += (rhs(i,j,k) - Ax*sig) / (s0*sig);
3367 }
3368 }
3369}
3370
3371AMREX_GPU_DEVICE AMREX_FORCE_INLINE
3372void mlndlap_gscolor_sten (int i, int j, int k, Array4<Real> const& sol,
3373 Array4<Real const> const& rhs,
3374 Array4<Real const> const& sten,
3375 Array4<int const> const& msk, int color) noexcept
3376{
3377 if (mlndlap_color(i,j,k) == color) {
3378 mlndlap_gauss_seidel_sten(i,j,k,sol,rhs,sten,msk);
3379 }
3380}
3381
3382}
3383#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_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
idir
int idir
Definition AMReX_HypreMLABecLap.cpp:1093
offset
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
fine
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
crse
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
AMREX_LOOP_3D
#define AMREX_LOOP_3D(bx, i, j, k, block)
Definition AMReX_Loop.nolint.H:4
amrex::HypreNodeLap::Int
HYPRE_Int Int
Definition AMReX_HypreNodeLap.H:36
amrex::Real
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
amrex::coarsen
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) coarsening ratio.
Definition AMReX_Box.H:1409
amrex
Definition AMReX_Amr.cpp:49
amrex::ubound
__host__ __device__ Dim3 ubound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:319
amrex::mlndlap_jacobi_ha
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_jacobi_ha(int i, int, int, Array4< Real > const &sol, Real Ax, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:92
amrex::mlndlap_set_connection
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_connection(int i, int j, int, Array4< Real > const &conn, Array4< Real const > const &intg, Array4< Real const > const &vol, Array4< EBCellFlag const > const &flag) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1689
amrex::mlndlap_adotx_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_aa(int i, int j, int k, Array4< Real const > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:66
amrex::mlndlap_mknewu_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_mknewu_c(int i, int, int, Array4< Real > const &u, Array4< Real const > const &p, Real sig, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:322
amrex::mlndlap_gscolor_ha
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_ha(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, int color) noexcept
Definition AMReX_MLNodeLap_1D_K.H:426
amrex::mlndlap_divu_fine_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_fine_contrib(int i, int j, int, Box const &fvbx, Box const &velbx, Array4< Real > const &rhs, Array4< Real const > const &vel, Array4< Real const > const &frhs, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, bool is_rz) noexcept
Definition AMReX_MLNodeLap_2D_K.H:867
amrex::surroundingNodes
__host__ __device__ BoxND< dim > surroundingNodes(const BoxND< dim > &b, int dir) noexcept
Return a BoxND with NODE based coordinates in direction dir that encloses BoxND b....
Definition AMReX_Box.H:1522
amrex::mlndlap_semi_avgdown_coeff
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_avgdown_coeff(int i, int j, int k, Array4< Real > const &crse, Array4< Real const > const &fine, int) noexcept
Definition AMReX_MLNodeLap_1D_K.H:30
amrex::mlndlap_gauss_seidel_aa
void mlndlap_gauss_seidel_aa(Box const &bx, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:193
amrex::mlndlap_interpadd_rap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_rap(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:408
amrex::mlndlap_avgdown_coeff_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_avgdown_coeff_y(int i, int j, int k, Array4< Real > const &crse, Array4< Real const > const &fine) noexcept
Definition AMReX_MLNodeLap_2D_K.H:36
amrex::mlndlap_gscolor_aa
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_aa(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, int color) noexcept
Definition AMReX_MLNodeLap_1D_K.H:448
amrex::mlndlap_normalize_ha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_normalize_ha(int i, int, int, Array4< Real > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:74
amrex::mlndlap_adotx_sten
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_sten(int, int, int, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:396
amrex::mlndlap_interpadd_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_c(int i, int, int, Array4< Real > const &fine, Array4< Real const > const &crse, Array4< int const > const &msk) noexcept
Definition AMReX_MLNodeLap_1D_K.H:233
amrex::mlndlap_set_stencil_s0
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_stencil_s0(int, int, int, Array4< Real > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:387
amrex::mlndlap_crse_resid
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_crse_resid(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &, Box const &, GpuArray< LinOpBCType, 3 > const &, GpuArray< LinOpBCType, 3 > const &, bool) noexcept
Definition AMReX_MLNodeLap_1D_K.H:349
amrex::mlndlap_gauss_seidel_with_line_solve_aa
void mlndlap_gauss_seidel_with_line_solve_aa(Box const &, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, GpuArray< Real, 3 > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:225
amrex::mlndlap_gscolor_c
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_c(int i, int j, int k, Array4< Real > const &sol, Array4< Real const > const &rhs, Real sig, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, int color) noexcept
Definition AMReX_MLNodeLap_1D_K.H:458
amrex::mlndlap_mknewu_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_mknewu_eb(int i, int j, int, Array4< Real > const &u, Array4< Real const > const &p, Array4< Real const > const &sig, Array4< Real const > const &vfrac, Array4< Real const > const &intg, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1818
amrex::mlndlap_Ax_fine_contrib_cs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_Ax_fine_contrib_cs(int i, int j, int, Box const &ndbx, Box const &ccbx, Array4< Real > const &f, Array4< Real const > const &res, Array4< Real const > const &rhs, Array4< Real const > const &phi, Real const sig, Array4< int const > const &msk, bool is_rz, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1157
amrex::mlndlap_res_cf_contrib_cs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_res_cf_contrib_cs(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Real, Array4< int const > const &, Array4< int const > const &, Array4< int const > const &, Array4< Real const > const &, GpuArray< Real, 3 > const &, Box const &, GpuArray< LinOpBCType, 3 > const &, GpuArray< LinOpBCType, 3 > const &, bool) noexcept
Definition AMReX_MLNodeLap_1D_K.H:369
amrex::Box
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:27
amrex::mlndlap_sum_Ax
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_sum_Ax(P const &pred, S const &sig, int i, int j, Real facx, Real facy, Array4< Real const > const &phi, bool is_rz) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1057
amrex::mlndlap_mknewu_eb_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_mknewu_eb_c(int i, int j, int, Array4< Real > const &u, Array4< Real const > const &p, Real sig, Array4< Real const > const &vfrac, Array4< Real const > const &intg, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1836
amrex::min
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:21
amrex::mlndlap_mknewu
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_mknewu(int i, int, int, Array4< Real > const &u, Array4< Real const > const &p, Array4< Real const > const &sig, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:312
amrex::mlndlap_rhcc_fine_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_rhcc_fine_contrib(int, int, int, Box const &, Array4< Real > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:332
amrex::mlndlap_semi_interpadd_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_semi_interpadd_aa(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, int) noexcept
Definition AMReX_MLNodeLap_1D_K.H:270
amrex::mlndlap_set_stencil
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_stencil(Box const &, Array4< Real > const &, Array4< Real const > const &, GpuArray< Real, 3 > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:381
amrex::neumann_scale
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real neumann_scale(int i, int j, Box const &nddom, GpuArray< LinOpBCType, 3 > const &bclo, GpuArray< LinOpBCType, 3 > const &bchi) noexcept
Definition AMReX_MLNodeLap_2D_K.H:930
amrex::aa_interp_line_y
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real aa_interp_line_y(Array4< Real const > const &crse, Array4< Real const > const &sig, int i, int j, int ic, int jc) noexcept
Definition AMReX_MLNodeLap_2D_K.H:582
amrex::mlndlap_set_integral_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_integral_eb(int i, int j, int, Array4< Real > const &intg, Array4< EBCellFlag const > const &flag, Array4< Real const > const &vol, Array4< Real const > const &ax, Array4< Real const > const &ay, Array4< Real const > const &bcen) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1894
amrex::mlndlap_zero_fine
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_zero_fine(int i, int, int, Array4< Real > const &phi, Array4< int const > const &msk, int fine_flag) noexcept
Definition AMReX_MLNodeLap_1D_K.H:8
amrex::mlndlap_Ax_fine_contrib_doit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_Ax_fine_contrib_doit(S const &sig, int i, int j, Box const &ndbx, Box const &ccbx, Array4< Real > const &f, Array4< Real const > const &res, Array4< Real const > const &rhs, Array4< Real const > const &phi, Array4< int const > const &msk, bool is_rz, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1099
amrex::mlndlap_set_stencil_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_stencil_eb(int i, int j, int, Array4< Real > const &sten, Array4< Real const > const &sig, Array4< Real const > const &conn, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1713
amrex::mlndlap_normalize_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_normalize_aa(int i, int j, int k, Array4< Real > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:84
amrex::ha_interp_face_xy
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real ha_interp_face_xy(Array4< Real const > const &crse, Array4< Real const > const &sigx, Array4< Real const > const &sigy, int i, int j, int ic, int jc) noexcept
Definition AMReX_MLNodeLap_2D_K.H:694
amrex::mlndlap_interpadd_aa
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_aa(int i, int, int, Array4< Real > const &fine, Array4< Real const > const &crse, Array4< Real const > const &sig, Array4< int const > const &msk) noexcept
Definition AMReX_MLNodeLap_1D_K.H:251
amrex::mlndlap_set_integral
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_integral(int i, int j, int, Array4< Real > const &intg) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1872
amrex::mlndlap_divu
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu(int i, int, int, Array4< Real > const &rhs, Array4< Real const > const &vel, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, Box const &, GpuArray< LinOpBCType, 3 > const &, GpuArray< LinOpBCType, 3 > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:284
amrex::mlndlap_Ax_fine_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_Ax_fine_contrib(int i, int j, int, Box const &ndbx, Box const &ccbx, Array4< Real > const &f, Array4< Real const > const &res, Array4< Real const > const &rhs, Array4< Real const > const &phi, Array4< Real const > const &sig, Array4< int const > const &msk, bool is_rz, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1143
amrex::mlndlap_gauss_seidel_sten
void mlndlap_gauss_seidel_sten(Box const &, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:401
amrex::mlndlap_rhcc_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_rhcc_eb(int i, int j, int, Array4< Real const > const &rhcc, Array4< Real const > const &vfrac, Array4< Real const > const &intg, Array4< int const > const &msk) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1854
amrex::mlndlap_adotx_ha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_ha(int i, int, int, Array4< Real const > const &x, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:52
amrex::mlndlap_jacobi_c
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_jacobi_c(int i, int, int, Array4< Real > const &sol, Real Ax, Array4< Real const > const &rhs, Real sig, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:141
amrex::mlndlap_set_surface_integral
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_surface_integral(int i, int j, int, Array4< Real > const &sintg) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1884
amrex::mlndlap_gscolor_sten
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_gscolor_sten(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, int) noexcept
Definition AMReX_MLNodeLap_1D_K.H:479
amrex::max
__host__ __device__ constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
amrex::mlndlap_gauss_seidel_ha
void mlndlap_gauss_seidel_ha(Box const &bx, Array4< Real > const &sol, Array4< Real const > const &rhs, Array4< Real const > const &sx, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:170
amrex::mlndlap_gauss_seidel_c
void mlndlap_gauss_seidel_c(Box const &bx, Array4< Real > const &sol, Array4< Real const > const &rhs, Real sig, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:203
amrex::mlndlap_adotx_c
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_c(int i, int, int, Array4< Real const > const &x, Real sigma, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:37
amrex::mlndlap_divu_cf_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_cf_contrib(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, Array4< int const > const &, Array4< int const > const &, GpuArray< Real, 3 > const &, Box const &, GpuArray< LinOpBCType, 3 > const &, GpuArray< LinOpBCType, 3 > const &, bool) noexcept
Definition AMReX_MLNodeLap_1D_K.H:338
amrex::mlndlap_restriction_rap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_restriction_rap(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:414
amrex::add_eb_flow_contribution
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void add_eb_flow_contribution(int i, int j, int, Array4< Real > const &rhs, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, Array4< Real const > const &bareaarr, Array4< Real const > const &sintg, Array4< Real const > const &eb_vel_dot_n) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1786
amrex::mlndlap_set_surface_integral_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_set_surface_integral_eb(int i, int j, int, Array4< Real > const &sintg, Array4< EBCellFlag const > const &flag, Array4< Real const > const &bcen, Array4< Real const > const &barea, Array4< Real const > const &bnorm) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1995
amrex::mlndlap_sum_Df
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_sum_Df(int ii, int jj, Real facx, Real facy, Array4< Real const > const &vel, Box const &velbx, bool is_rz) noexcept
Definition AMReX_MLNodeLap_2D_K.H:843
amrex::LoopConcurrentOnCpu
void LoopConcurrentOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:388
amrex::aa_interp_face_xy
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real aa_interp_face_xy(Array4< Real const > const &crse, Array4< Real const > const &sig, int i, int j, int ic, int jc) noexcept
Definition AMReX_MLNodeLap_2D_K.H:591
amrex::mlndlap_color
AMREX_GPU_DEVICE AMREX_FORCE_INLINE int mlndlap_color(int i, int, int)
Definition AMReX_MLNodeLap_1D_K.H:420
amrex::Abort
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:230
amrex::mlndlap_jacobi_aa
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlndlap_jacobi_aa(int i, int j, int k, Array4< Real > const &sol, Real Ax, Array4< Real const > const &rhs, Array4< Real const > const &sig, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv) noexcept
Definition AMReX_MLNodeLap_1D_K.H:125
amrex::mlndlap_avgdown_coeff_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_avgdown_coeff_x(int i, int, int, Array4< Real > const &crse, Array4< Real const > const &fine) noexcept
Definition AMReX_MLNodeLap_1D_K.H:21
amrex::LoopOnCpu
void LoopOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:365
amrex::mlndlap_interpadd_ha
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_interpadd_ha(int i, int j, int k, Array4< Real > const &fine, Array4< Real const > const &crse, Array4< Real const > const &sig, Array4< int const > const &msk) noexcept
Definition AMReX_MLNodeLap_1D_K.H:276
amrex::aa_interp_line_x
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real aa_interp_line_x(Array4< Real const > const &crse, Array4< Real const > const &sig, int i, int j, int ic, int jc) noexcept
Definition AMReX_MLNodeLap_2D_K.H:573
amrex::mlndlap_stencil_rap
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_stencil_rap(int, int, int, Array4< Real > const &, Array4< Real const > const &) noexcept
Definition AMReX_MLNodeLap_1D_K.H:391
amrex::mlndlap_adotx_sten_doit
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_adotx_sten_doit(int i, int j, int k, Array4< Real const > const &x, Array4< Real const > const &sten) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1507
amrex::tridiagonal_solve
AMREX_FORCE_INLINE void tridiagonal_solve(Array1D< T, 0, 31 > &a_ls, Array1D< T, 0, 31 > &b_ls, Array1D< T, 0, 31 > &c_ls, Array1D< T, 0, 31 > &r_ls, Array1D< T, 0, 31 > &u_ls, Array1D< T, 0, 31 > &gam, int ilen) noexcept
Definition AMReX_MLABecLap_3D_K.H:432
amrex::mlndlap_rhcc
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real mlndlap_rhcc(int i, int, int, Array4< Real const > const &rhcc, Array4< int const > const &msk) noexcept
Definition AMReX_MLNodeLap_1D_K.H:299
amrex::scale
__host__ __device__ IntVectND< dim > scale(const IntVectND< dim > &p, int s) noexcept
Returns a IntVectND obtained by multiplying each of the components of this IntVectND by s.
Definition AMReX_IntVect.H:1016
amrex::LoopConcurrent
__host__ __device__ void LoopConcurrent(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:152
amrex::mlndlap_divu_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_divu_eb(int i, int j, int, Array4< Real > const &rhs, Array4< Real const > const &vel, Array4< Real const > const &vfrac, Array4< Real const > const &intg, Array4< int const > const &msk, GpuArray< Real, 3 > const &dxinv, Box const &nodal_domain, GpuArray< LinOpBCType, 3 > const &bclo, GpuArray< LinOpBCType, 3 > const &bchi) noexcept
Definition AMReX_MLNodeLap_2D_K.H:1729
amrex::lbound
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Definition AMReX_Array4.H:312
amrex::mlndlap_res_cf_contrib
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mlndlap_res_cf_contrib(int, int, int, Array4< Real > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< Real const > const &, Array4< int const > const &, Array4< int const > const &, Array4< int const > const &, Array4< Real const > const &, GpuArray< Real, 3 > const &, Box const &, GpuArray< LinOpBCType, 3 > const &, GpuArray< LinOpBCType, 3 > const &, bool) noexcept
Definition AMReX_MLNodeLap_1D_K.H:357
amrex::Array1D
Definition AMReX_Array.H:199
amrex::Array4
Definition AMReX_Array4.H:61
amrex::GpuArray
Fixed-size array that can be used on GPU.
Definition AMReX_Array.H:40