Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLCurlCurl_K.H
Go to the documentation of this file.
1#ifndef AMREX_ML_CURL_CURL_K_H_
2#define AMREX_ML_CURL_CURL_K_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_Array4.H>
6#include <AMReX_LUSolver.H>
7#include <AMReX_PCGSolver.H>
8
9namespace amrex {
10
11/* Index types
12 * E_x : (0,1,1)
13 * E_y : (1,0,1)
14 * E_z : (1,1,0)
15 * (curl E)_x : (1,0,0)
16 * (curl E)_y : (0,1,0)
17 * (curl E)_z : (0,0,1)
18 */
19
20/*
21 Notes for gs4:
22
23 Interior nodes:
24
25 For v = [Ex(i-1,j,k), Ex(i,j,k), Ey(i,j-1,k), Ey(i,j,k), Ez(i,j,k-1), Ez(i,j,k)]^T,
26 we have A*v = b, where
27
28 A00 = alpha*(dyy+dzz)*2 + beta
29 A01 = 0
30 A02 = -alpha*dxy
31 A03 = alpha*dxy
32 A04 = -alpha*dxz
33 A05 = alpha*dxz
34
35 A10 = 0
36 A11 = alpha*(dyy+dzz)*2 + beta
37 A12 = alpha*dxy
38 A13 = -alpha*dxy
39 A14 = alpha*dxz
40 A15 = -alpha*dxz
41
42 A20 = -alpha*dxy
43 A21 = alpha*dxy
44 A22 = alpha*(dxx+dzz)*2 + beta
45 A23 = 0
46 A24 = -alpha*dyz
47 A25 = alpha*dyz
48
49 A30 = alpha*dxy
50 A31 = -alpha*dxy
51 A32 = 0
52 A33 = alpha*(dxx+dzz)*2 + beta
53 A34 = alpha*dyz
54 A35 = -alpha*dyz
55
56 A40 = -alpha*dxz
57 A41 = alpha*dxz
58 A42 = -alpha*dyz
59 A43 = alpha*dyz
60 A44 = alpha*(dxx+dyy)*2 + beta
61 A45 = 0
62
63 A50 = alpha*dxz
64 A51 = -alpha*dxz
65 A52 = alpha*dyz
66 A53 = -alpha*dyz
67 A54 = 0
68 A55 = alpha*(dxx+dyy)*2 + beta
69
70 b0 = rhsx(i-1,j,k) - (alpha*ccex), where
71 ccex = - dyy * (ex(i-1,j-1,k ) +
72 ex(i-1,j+1,k ))
73 - dzz * (ex(i-1,j ,k+1) +
74 ex(i-1,j ,k-1))
75 + dxy * (ey(i-1,j-1,k )
76 - ey(i-1,j ,k ))
77 + dxz * (ez(i-1,j ,k-1)
78 - ez(i-1,j ,k ))
79 b1 = rhsx(i,j,k) - (alpha*ccex), where
80 ccex = - dyy * ( ex(i ,j-1,k ) +
81 ex(i ,j+1,k ))
82 - dzz * ( ex(i ,j ,k+1) +
83 ex(i ,j ,k-1))
84 + dxy * (-ey(i+1,j-1,k )
85 + ey(i+1,j ,k ))
86 + dxz * (-ez(i+1,j ,k-1)
87 + ez(i+1,j ,k ));
88 b2 = rhsy(i,j-1,k) - alpha*ccey, where
89 ccey = - dxx * (ey(i-1,j-1,k ) +
90 ey(i+1,j-1,k ))
91 - dzz * (ey(i ,j-1,k-1) +
92 ey(i ,j-1,k+1))
93 + dxy * (ex(i-1,j-1,k )
94 - ex(i ,j-1,k ))
95 + dyz * (ez(i ,j-1,k-1)
96 - ez(i ,j-1,k ))
97 b3 = rhsy(i,j,k) - alpha*ccey, where
98 ccey = - dxx * ( ey(i-1,j ,k ) +
99 ey(i+1,j ,k ))
100 - dzz * ( ey(i ,j ,k-1) +
101 ey(i ,j ,k+1))
102 + dxy * (-ex(i-1,j+1,k )
103 + ex(i ,j+1,k ))
104 + dyz * (-ez(i ,j+1,k-1)
105 + ez(i ,j+1,k ));
106 b4 = rhsz(i,j,k-1) - alpha*ccez, where
107 ccez = - dxx * (ez(i-1,j ,k-1) +
108 ez(i+1,j ,k-1))
109 - dyy * (ez(i ,j-1,k-1) +
110 ez(i ,j+1,k-1))
111 + dxz * (ex(i-1,j ,k-1)
112 - ex(i ,j ,k-1))
113 + dyz * (ey(i ,j-1,k-1)
114 - ey(i ,j ,k-1))
115 b5 = rhsz(i,j,k) - alpha*ccez, where
116 ccez = - dxx * ( ez(i-1,j ,k ) +
117 ez(i+1,j ,k ))
118 - dyy * ( ez(i ,j-1,k ) +
119 ez(i ,j+1,k ))
120 + dxz * (-ex(i-1,j ,k+1)
121 + ex(i ,j ,k+1))
122 + dyz * (-ey(i ,j-1,k+1)
123 + ey(i ,j ,k+1));
124
125 dxx = 1/(dx*dx)
126 dyy = 1/(dy*dy)
127 dzz = 1/(dz*dz)
128 dxy = 1/(dx*dy)
129 dxz = 1/(dx*dz)
130 dyz = 1/(dy*dz)
131
132 For Dirichlet boundary nodes, we don't do anything.
133
134 For symmetric boundary nodes, we treat it as interior nodes because the
135 rhs outside the domain has been filled properly.
136
137 In 2D,
138
139 For v = [Ex(i-1,j,k), Ex(i,j,k), Ey(i,j-1,k), Ey(i,j,k)]^T,
140 we have A*v = b, where
141
142 A00 = alpha*dyy*2 + beta
143 A01 = 0
144 A02 = -alpha*dxy
145 A03 = alpha*dxy
146
147 A10 = 0
148 A11 = alpha*dyy*2 + beta
149 A12 = alpha*dxy
150 A13 = -alpha*dxy
151
152 A20 = -alpha*dxy
153 A21 = alpha*dxy
154 A22 = alpha*dxx*2 + beta
155 A23 = 0
156
157 A30 = alpha*dxy
158 A31 = -alpha*dxy
159 A32 = 0
160 A33 = alpha*dxx*2 + beta
161
162 b0 = rhsx(i-1,j,k) - (alpha*ccex), where
163 ccex = - dyy * (ex(i-1,j-1,k ) +
164 ex(i-1,j+1,k ))
165 + dxy * (ey(i-1,j-1,k )
166 - ey(i-1,j ,k ))
167 b1 = rhsx(i,j,k) - (alpha*ccex), where
168 ccex = - dyy * ( ex(i ,j-1,k ) +
169 ex(i ,j+1,k ))
170 + dxy * (-ey(i+1,j-1,k )
171 + ey(i+1,j ,k ))
172 b2 = rhsy(i,j-1,k) - alpha*ccey, where
173 ccey = - dxx * (ey(i-1,j-1,k ) +
174 ey(i+1,j-1,k ))
175 + dxy * (ex(i-1,j-1,k )
176 - ex(i ,j-1,k ))
177 b3 = rhsy(i,j,k) - alpha*ccey, where
178 ccey = - dxx * ( ey(i-1,j ,k ) +
179 ey(i+1,j ,k ))
180 + dxy * (-ex(i-1,j+1,k )
181 + ex(i ,j+1,k ))
182*/
183
184struct CurlCurlDirichletInfo
185{
186 IntVect dirichlet_lo;
187 IntVect dirichlet_hi;
188
189 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
190 bool is_dirichlet_node (int i, int j, int k) const
191 {
192#if (AMREX_SPACEDIM == 1)
193 amrex::ignore_unused(j,k);
194 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0]);
195#elif (AMREX_SPACEDIM == 2)
196 amrex::ignore_unused(k);
197 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
198 || (j == dirichlet_lo[1]) || (j == dirichlet_hi[1]);
199#else
200 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
201 || (j == dirichlet_lo[1]) || (j == dirichlet_hi[1])
202 || (k == dirichlet_lo[2]) || (k == dirichlet_hi[2]);
203#endif
204 }
205
206 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
207 bool is_dirichlet_x_edge (int, int j, int k) const // NOLINT(readability-convert-member-functions-to-static)
208 {
209#if (AMREX_SPACEDIM == 1)
210 amrex::ignore_unused(j, k);
211 return false;
212#elif (AMREX_SPACEDIM == 2)
213 amrex::ignore_unused(k);
214 return (j == dirichlet_lo[1]) || (j == dirichlet_hi[1]);
215#else
216 return (j == dirichlet_lo[1]) || (j == dirichlet_hi[1])
217 || (k == dirichlet_lo[2]) || (k == dirichlet_hi[2]);
218#endif
219 }
220
221 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
222 bool is_dirichlet_y_edge (int i, int, int k) const
223 {
224#if (AMREX_SPACEDIM < 3)
225 amrex::ignore_unused(k);
226 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0]);
227#else
228 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
229 || (k == dirichlet_lo[2]) || (k == dirichlet_hi[2]);
230#endif
231 }
232
233 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
234 bool is_dirichlet_z_edge (int i, int j, int) const
235 {
236#if (AMREX_SPACEDIM == 1)
237 amrex::ignore_unused(j);
238 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0]);
239#else
240 return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
241 || (j == dirichlet_lo[1]) || (j == dirichlet_hi[1]);
242#endif
243 }
244
245 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
246 bool is_dirichlet_edge (int dim, int i, int j, int k) const
247 {
248 if (dim == 0) {
249 return is_dirichlet_x_edge(i,j,k);
250 } else if (dim == 1) {
251 return is_dirichlet_y_edge(i,j,k);
252 } else {
253 return is_dirichlet_z_edge(i,j,k);
254 }
255 }
256};
257
258struct CurlCurlSymmetryInfo
259{
260 IntVect symmetry_lo;
261 IntVect symmetry_hi;
262
263 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
264 bool xlo_is_symmetric (int i) const
265 {
266 return i == symmetry_lo[0];
267 }
268
269 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
270 bool xhi_is_symmetric (int i) const
271 {
272 return i == symmetry_hi[0];
273 }
274
275#if (AMREX_SPACEDIM > 1)
276 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
277 bool ylo_is_symmetric (int j) const
278 {
279 return j == symmetry_lo[1];
280 }
281
282 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
283 bool yhi_is_symmetric (int j) const
284 {
285 return j == symmetry_hi[1];
286 }
287#endif
288
289#if (AMREX_SPACEDIM == 3)
290 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
291 bool zlo_is_symmetric (int k) const
292 {
293 return k == symmetry_lo[2];
294 }
295
296 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
297 bool zhi_is_symmetric (int k) const
298 {
299 return k == symmetry_hi[2];
300 }
301#endif
302
303 [[nodiscard]] bool is_symmetric (int dir, int side, int idx) const
304 {
305#if (AMREX_SPACEDIM == 1)
306 amrex::ignore_unused(dir);
307 return (side == 0) ? xlo_is_symmetric(idx) : xhi_is_symmetric(idx);
308#elif (AMREX_SPACEDIM == 2)
309 if (dir == 0) {
310 return (side == 0) ? xlo_is_symmetric(idx) : xhi_is_symmetric(idx);
311 } else {
312 return (side == 0) ? ylo_is_symmetric(idx) : yhi_is_symmetric(idx);
313 }
314#else
315 if (dir == 0) {
316 return (side == 0) ? xlo_is_symmetric(idx) : xhi_is_symmetric(idx);
317 } else if (dir == 1) {
318 return (side == 0) ? ylo_is_symmetric(idx) : yhi_is_symmetric(idx);
319 } else {
320 return (side == 0) ? zlo_is_symmetric(idx) : zhi_is_symmetric(idx);
321 }
322#endif
323 }
324};
325
326enum struct CurlCurlStateType { x, b, r }; // x, b & r=b-Ax
327
328AMREX_GPU_DEVICE AMREX_FORCE_INLINE
329void mlcurlcurl_adotx_x (int i, int j, int k, Array4<Real> const& Ax,
330 Array4<Real const> const& ex,
331 Array4<Real const> const& ey,
332 Array4<Real const> const& ez,
333 Real beta, GpuArray<Real,AMREX_SPACEDIM> const& adxinv)
334{
335#if (AMREX_SPACEDIM == 1)
336 amrex::ignore_unused(ey,ez,adxinv);
337 Real ccex = 0;
338#elif (AMREX_SPACEDIM == 2)
339 amrex::ignore_unused(ez);
340 Real dyy = adxinv[1] * adxinv[1];
341 Real dxy = adxinv[0] * adxinv[1];
342 Real ccex = ex(i ,j ,k ) * dyy * Real(2.0)
343 - dyy * (ex(i ,j-1,k ) +
344 ex(i ,j+1,k ))
345 + dxy * (ey(i ,j-1,k )
346 - ey(i ,j ,k )
347 - ey(i+1,j-1,k )
348 + ey(i+1,j ,k ));
349#elif (AMREX_SPACEDIM == 3)
350 Real dyy = adxinv[1] * adxinv[1];
351 Real dzz = adxinv[2] * adxinv[2];
352 Real dxy = adxinv[0] * adxinv[1];
353 Real dxz = adxinv[0] * adxinv[2];
354 Real ccex = ex(i ,j ,k ) * (dyy+dzz)*Real(2.0)
355 - dyy * (ex(i ,j-1,k ) +
356 ex(i ,j+1,k ))
357 - dzz * (ex(i ,j ,k+1) +
358 ex(i ,j ,k-1))
359 + dxy * (ey(i ,j-1,k )
360 - ey(i ,j ,k )
361 - ey(i+1,j-1,k )
362 + ey(i+1,j ,k ))
363 + dxz * (ez(i ,j ,k-1)
364 - ez(i ,j ,k )
365 - ez(i+1,j ,k-1)
366 + ez(i+1,j ,k ));
367#endif
368 Ax(i,j,k) = ccex + beta * ex(i,j,k);
369}
370
371AMREX_GPU_DEVICE AMREX_FORCE_INLINE
372void mlcurlcurl_adotx_y (int i, int j, int k, Array4<Real> const& Ay,
373 Array4<Real const> const& ex,
374 Array4<Real const> const& ey,
375 Array4<Real const> const& ez,
376 Real beta, GpuArray<Real,AMREX_SPACEDIM> const& adxinv)
377{
378#if (AMREX_SPACEDIM == 1)
379 amrex::ignore_unused(ex,ez);
380 Real dxx = adxinv[0] * adxinv[0];
381 Real ccey = ey(i ,j ,k ) * dxx * Real(2.0)
382 - dxx * (ey(i-1,j ,k ) +
383 ey(i+1,j ,k ));
384#elif (AMREX_SPACEDIM == 2)
385 amrex::ignore_unused(ez);
386 Real dxx = adxinv[0] * adxinv[0];
387 Real dxy = adxinv[0] * adxinv[1];
388 Real ccey = ey(i ,j ,k ) * dxx * Real(2.0)
389 - dxx * (ey(i-1,j ,k ) +
390 ey(i+1,j ,k ))
391 + dxy * (ex(i-1,j ,k )
392 - ex(i ,j ,k )
393 - ex(i-1,j+1,k )
394 + ex(i ,j+1,k ));
395#elif (AMREX_SPACEDIM == 3)
396 Real dxx = adxinv[0] * adxinv[0];
397 Real dzz = adxinv[2] * adxinv[2];
398 Real dxy = adxinv[0] * adxinv[1];
399 Real dyz = adxinv[1] * adxinv[2];
400 Real ccey = ey(i ,j ,k ) * (dxx+dzz)*Real(2.0)
401 - dxx * (ey(i-1,j ,k ) +
402 ey(i+1,j ,k ))
403 - dzz * (ey(i ,j ,k-1) +
404 ey(i ,j ,k+1))
405 + dxy * (ex(i-1,j ,k )
406 - ex(i ,j ,k )
407 - ex(i-1,j+1,k )
408 + ex(i ,j+1,k ))
409 + dyz * (ez(i ,j ,k-1)
410 - ez(i ,j ,k )
411 - ez(i ,j+1,k-1)
412 + ez(i ,j+1,k ));
413#endif
414 Ay(i,j,k) = ccey + beta * ey(i,j,k);
415}
416
417AMREX_GPU_DEVICE AMREX_FORCE_INLINE
418void mlcurlcurl_adotx_z (int i, int j, int k, Array4<Real> const& Az,
419 Array4<Real const> const& ex,
420 Array4<Real const> const& ey,
421 Array4<Real const> const& ez,
422 Real beta, GpuArray<Real,AMREX_SPACEDIM> const& adxinv)
423{
424#if (AMREX_SPACEDIM == 1)
425 amrex::ignore_unused(ex,ey);
426 Real dxx = adxinv[0] * adxinv[0];
427 Real ccez = ez(i ,j ,k ) * dxx*Real(2.0)
428 - dxx * (ez(i-1,j ,k ) +
429 ez(i+1,j ,k ));
430#elif (AMREX_SPACEDIM == 2)
431 amrex::ignore_unused(ex,ey);
432 Real dxx = adxinv[0] * adxinv[0];
433 Real dyy = adxinv[1] * adxinv[1];
434 Real ccez = ez(i ,j ,k ) * (dxx+dyy)*Real(2.0)
435 - dxx * (ez(i-1,j ,k ) +
436 ez(i+1,j ,k ))
437 - dyy * (ez(i ,j-1,k ) +
438 ez(i ,j+1,k ));
439#elif (AMREX_SPACEDIM == 3)
440 Real dxx = adxinv[0] * adxinv[0];
441 Real dyy = adxinv[1] * adxinv[1];
442 Real dxz = adxinv[0] * adxinv[2];
443 Real dyz = adxinv[1] * adxinv[2];
444 Real ccez = ez(i ,j ,k ) * (dxx+dyy)*Real(2.0)
445 - dxx * (ez(i-1,j ,k ) +
446 ez(i+1,j ,k ))
447 - dyy * (ez(i ,j-1,k ) +
448 ez(i ,j+1,k ))
449 + dxz * (ex(i-1,j ,k )
450 - ex(i ,j ,k )
451 - ex(i-1,j ,k+1)
452 + ex(i ,j ,k+1))
453 + dyz * (ey(i ,j-1,k )
454 - ey(i ,j ,k )
455 - ey(i ,j-1,k+1)
456 + ey(i ,j ,k+1));
457#endif
458 Az(i,j,k) = ccez + beta * ez(i,j,k);
459}
460
461#if (AMREX_SPACEDIM == 1)
462AMREX_GPU_DEVICE AMREX_FORCE_INLINE
463void mlcurlcurl_1D (int i, int j, int k,
464 Array4<Real> const& ex,
465 Array4<Real> const& ey,
466 Array4<Real> const& ez,
467 Array4<Real const> const& rhsx,
468 Array4<Real const> const& rhsy,
469 Array4<Real const> const& rhsz,
470 Real beta,
471 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
472 int color,
473 CurlCurlDirichletInfo const& dinfo)
474{
475 if (dinfo.is_dirichlet_node(i,j,k)) {return; }
476
477 Real dxx = adxinv[0] * adxinv[0];
478
479 int my_color = i%2;
480
481 if (my_color == color)
482 {
483 ex(i,j,k) = rhsx(i,j,k) / beta;
484
485 Real gamma_y = dxx * Real(2.0) + beta;
486 Real ccey = - dxx * (ey(i-1,j ,k ) +
487 ey(i+1,j ,k ) );
488 Real res_y = rhsy(i,j,k) - ( gamma_y * ey(i,j,k) + ccey );
489 ey(i,j,k) += res_y/gamma_y;
490
491 Real gamma_z = dxx * Real(2.0) + beta;
492 Real ccez = -dxx * (ez(i-1,j ,k ) +
493 ez(i+1,j ,k ) );
494 Real res_z = rhsz(i,j,k) - ( gamma_z * ez(i,j,k) + ccez );
495 ez(i,j,k) += res_z/gamma_z;
496 }
497}
498
499
500AMREX_GPU_DEVICE AMREX_FORCE_INLINE
501void mlcurlcurl_1D (int i, int j, int k,
502 Array4<Real> const& ex,
503 Array4<Real> const& ey,
504 Array4<Real> const& ez,
505 Array4<Real const> const& rhsx,
506 Array4<Real const> const& rhsy,
507 Array4<Real const> const& rhsz,
508 Array4<Real const> const& betax,
509 Array4<Real const> const& betay,
510 Array4<Real const> const& betaz,
511 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
512 int color,
513 CurlCurlDirichletInfo const& dinfo)
514{
515 if (dinfo.is_dirichlet_node(i,j,k)) {return; }
516 Real dxx = adxinv[0] * adxinv[0];
517
518 int my_color = i%2;
519
520 if (my_color == color)
521 {
522 ex(i,j,k) = rhsx(i,j,k) / betax(i,j,k);
523
524 Real gamma_y = dxx * Real(2.0) + betay(i,j,k);
525 Real ccey = - dxx * (ey(i-1,j ,k ) +
526 ey(i+1,j ,k ) );
527 Real res_y = rhsy(i,j,k) - ( gamma_y * ey(i,j,k) + ccey );
528 ey(i,j,k) += res_y/gamma_y;
529
530 Real gamma_z = dxx * Real(2.0) + betaz(i,j,k);
531 Real ccez = -dxx * (ez(i-1,j ,k ) +
532 ez(i+1,j ,k ) );
533 Real res_z = rhsz(i,j,k) - ( gamma_z * ez(i,j,k) + ccez );
534 ez(i,j,k) += res_z/gamma_z;
535 }
536}
537
538
539#endif
540
541#if (AMREX_SPACEDIM > 1)
542
543AMREX_GPU_DEVICE AMREX_FORCE_INLINE
544void mlcurlcurl_gs4_lu (int i, int j, int k,
545 Array4<Real> const& ex,
546 Array4<Real> const& ey,
547 Array4<Real> const& ez,
548 Array4<Real const> const& rhsx,
549 Array4<Real const> const& rhsy,
550 Array4<Real const> const& rhsz,
551#if (AMREX_SPACEDIM == 2)
552 Real beta,
553#endif
554 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
555 int color, LUSolver<AMREX_SPACEDIM*2,Real> const& lusolver,
556 CurlCurlDirichletInfo const& dinfo,
557 CurlCurlSymmetryInfo const& sinfo)
558{
559 if (dinfo.is_dirichlet_node(i,j,k)) { return; }
560
561 int my_color = i%2 + 2*(j%2);
562 if (k%2 != 0) {
563 my_color = 3 - my_color;
564 }
565
566#if (AMREX_SPACEDIM == 2)
567
568 Real dxx = adxinv[0] * adxinv[0];
569 Real dyy = adxinv[1] * adxinv[1];
570
571 if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) ||
572 ((my_color == 1 || my_color == 2) && (color == 1 || color == 2)))
573 {
574 Real gamma = (dxx+dyy)*Real(2.0) + beta;
575 Real ccez = - dxx * (ez(i-1,j ,k ) +
576 ez(i+1,j ,k ))
577 - dyy * (ez(i ,j-1,k ) +
578 ez(i ,j+1,k ));
579 Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez);
580 constexpr Real omega = Real(1.15);
581 ez(i,j,k) += omega/gamma * res;
582 }
583
584 if (my_color != color) { return; }
585
586 Real dxy = adxinv[0]*adxinv[1];
587
588 GpuArray<Real,6> b
589 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
590 ex(i-1,j+1,k ))
591 + dxy * ( ey(i-1,j-1,k )
592 -ey(i-1,j ,k ))),
593 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
594 ex(i ,j+1,k ))
595 + dxy * (-ey(i+1,j-1,k )
596 +ey(i+1,j ,k ))),
597 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
598 ey(i+1,j-1,k ))
599 + dxy * ( ex(i-1,j-1,k )
600 -ex(i ,j-1,k ))),
601 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
602 ey(i+1,j ,k ))
603 + dxy * (-ex(i-1,j+1,k )
604 +ex(i ,j+1,k )))};
605
606 if (sinfo.xlo_is_symmetric(i)) {
607 b[0] = -b[1];
608 } else if (sinfo.xhi_is_symmetric(i)) {
609 b[1] = -b[0];
610 }
611
612 if (sinfo.ylo_is_symmetric(j)) {
613 b[2] = -b[3];
614 } else if (sinfo.yhi_is_symmetric(j)) {
615 b[3] = -b[2];
616 }
617
618 GpuArray<Real,4> x;
619 lusolver(x.data(), b.data());
620 ex(i-1,j ,k ) = x[0];
621 ex(i ,j ,k ) = x[1];
622 ey(i ,j-1,k ) = x[2];
623 ey(i ,j ,k ) = x[3];
624
625#elif (AMREX_SPACEDIM == 3)
626
627 if (my_color != color) { return; }
628
629 Real dxx = adxinv[0]*adxinv[0];
630 Real dyy = adxinv[1]*adxinv[1];
631 Real dzz = adxinv[2]*adxinv[2];
632 Real dxy = adxinv[0]*adxinv[1];
633 Real dxz = adxinv[0]*adxinv[2];
634 Real dyz = adxinv[1]*adxinv[2];
635
636 GpuArray<Real,6> b
637 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
638 ex(i-1,j+1,k ))
639 - dzz * ( ex(i-1,j ,k+1) +
640 ex(i-1,j ,k-1))
641 + dxy * ( ey(i-1,j-1,k )
642 -ey(i-1,j ,k ))
643 + dxz * ( ez(i-1,j ,k-1)
644 -ez(i-1,j ,k ))),
645 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
646 ex(i ,j+1,k ))
647 - dzz * ( ex(i ,j ,k+1) +
648 ex(i ,j ,k-1))
649 + dxy * (-ey(i+1,j-1,k )
650 +ey(i+1,j ,k ))
651 + dxz * (-ez(i+1,j ,k-1)
652 +ez(i+1,j ,k ))),
653 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
654 ey(i+1,j-1,k ))
655 - dzz * ( ey(i ,j-1,k-1) +
656 ey(i ,j-1,k+1))
657 + dxy * ( ex(i-1,j-1,k )
658 -ex(i ,j-1,k ))
659 + dyz * ( ez(i ,j-1,k-1)
660 -ez(i ,j-1,k ))),
661 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
662 ey(i+1,j ,k ))
663 - dzz * ( ey(i ,j ,k-1) +
664 ey(i ,j ,k+1))
665 + dxy * (-ex(i-1,j+1,k )
666 +ex(i ,j+1,k ))
667 + dyz * (-ez(i ,j+1,k-1)
668 +ez(i ,j+1,k ))),
669 rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) +
670 ez(i+1,j ,k-1))
671 - dyy * ( ez(i ,j-1,k-1) +
672 ez(i ,j+1,k-1))
673 + dxz * ( ex(i-1,j ,k-1)
674 -ex(i ,j ,k-1))
675 + dyz * ( ey(i ,j-1,k-1)
676 -ey(i ,j ,k-1))),
677 rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) +
678 ez(i+1,j ,k ))
679 - dyy * ( ez(i ,j-1,k ) +
680 ez(i ,j+1,k ))
681 + dxz * (-ex(i-1,j ,k+1)
682 +ex(i ,j ,k+1))
683 + dyz * (-ey(i ,j-1,k+1)
684 +ey(i ,j ,k+1)))};
685
686 if (sinfo.xlo_is_symmetric(i)) {
687 b[0] = -b[1];
688 } else if (sinfo.xhi_is_symmetric(i)) {
689 b[1] = -b[0];
690 }
691
692 if (sinfo.ylo_is_symmetric(j)) {
693 b[2] = -b[3];
694 } else if (sinfo.yhi_is_symmetric(j)) {
695 b[3] = -b[2];
696 }
697
698 if (sinfo.zlo_is_symmetric(k)) {
699 b[4] = -b[5];
700 } else if (sinfo.zhi_is_symmetric(k)) {
701 b[5] = -b[4];
702 }
703
704 GpuArray<Real,6> x;
705 lusolver(x.data(), b.data());
706 ex(i-1,j ,k ) = x[0];
707 ex(i ,j ,k ) = x[1];
708 ey(i ,j-1,k ) = x[2];
709 ey(i ,j ,k ) = x[3];
710 ez(i ,j ,k-1) = x[4];
711 ez(i ,j ,k ) = x[5];
712#endif
713}
714
715template <bool PCG>
716AMREX_GPU_DEVICE AMREX_FORCE_INLINE
717void mlcurlcurl_gs4 (int i, int j, int k,
718 Array4<Real> const& ex,
719 Array4<Real> const& ey,
720 Array4<Real> const& ez,
721 Array4<Real const> const& rhsx,
722 Array4<Real const> const& rhsy,
723 Array4<Real const> const& rhsz,
724 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
725 int color,
726 Array4<Real const> const& betax,
727 Array4<Real const> const& betay,
728 Array4<Real const> const& betaz,
729 CurlCurlDirichletInfo const& dinfo,
730 CurlCurlSymmetryInfo const& sinfo)
731{
732 if (dinfo.is_dirichlet_node(i,j,k)) { return; }
733
734 int my_color = i%2 + 2*(j%2);
735 if (k%2 != 0) {
736 my_color = 3 - my_color;
737 }
738
739#if (AMREX_SPACEDIM == 2)
740
741 Real dxx = adxinv[0] * adxinv[0];
742 Real dyy = adxinv[1] * adxinv[1];
743
744 if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) ||
745 ((my_color == 1 || my_color == 2) && (color == 1 || color == 2)))
746 {
747 Real gamma = (dxx+dyy)*Real(2.0) + betaz(i,j,k);
748 Real ccez = - dxx * (ez(i-1,j ,k ) +
749 ez(i+1,j ,k ))
750 - dyy * (ez(i ,j-1,k ) +
751 ez(i ,j+1,k ));
752 Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez);
753 constexpr Real omega = Real(1.15);
754 ez(i,j,k) += omega/gamma * res;
755 }
756
757 if (my_color != color) { return; }
758
759 Real dxy = adxinv[0]*adxinv[1];
760
761 GpuArray<Real,6> b
762 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
763 ex(i-1,j+1,k ))
764 + dxy * ( ey(i-1,j-1,k )
765 -ey(i-1,j ,k ))),
766 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
767 ex(i ,j+1,k ))
768 + dxy * (-ey(i+1,j-1,k )
769 +ey(i+1,j ,k ))),
770 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
771 ey(i+1,j-1,k ))
772 + dxy * ( ex(i-1,j-1,k )
773 -ex(i ,j-1,k ))),
774 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
775 ey(i+1,j ,k ))
776 + dxy * (-ex(i-1,j+1,k )
777 +ex(i ,j+1,k )))};
778
779 GpuArray<Real,4> beta;
780
781 if (sinfo.xlo_is_symmetric(i)) {
782 b[0] = -b[1];
783 beta[0] = beta[1] = betax(i,j,k);
784 } else if (sinfo.xhi_is_symmetric(i)) {
785 b[1] = -b[0];
786 beta[0] = beta[1] = betax(i-1,j,k);
787 } else {
788 beta[0] = betax(i-1,j,k);
789 beta[1] = betax(i ,j,k);
790 }
791
792 if (sinfo.ylo_is_symmetric(j)) {
793 b[2] = -b[3];
794 beta[2] = beta[3] = betay(i,j,k);
795 } else if (sinfo.yhi_is_symmetric(j)) {
796 b[3] = -b[2];
797 beta[2] = beta[3] = betay(i,j-1,k);
798 } else {
799 beta[2] = betay(i,j-1,k);
800 beta[3] = betay(i,j ,k);
801 }
802
803 if constexpr (PCG) {
804 Real diagInv[4] = {Real(1.0) / (dyy*Real(2.0) + beta[0]),
805 Real(1.0) / (dyy*Real(2.0) + beta[1]),
806 Real(1.0) / (dxx*Real(2.0) + beta[2]),
807 Real(1.0) / (dxx*Real(2.0) + beta[3])};
808 auto precond = [&] (Real * AMREX_RESTRICT z,
809 Real const* AMREX_RESTRICT r)
810 {
811 for (int m = 0; m < 4; ++m) { z[m] = r[m] * diagInv[m]; }
812 };
813 auto mat = [&] (Real * AMREX_RESTRICT Av,
814 Real const* AMREX_RESTRICT v)
815 {
816 Av[0] = (dyy*Real(2.0) + beta[0]) * v[0] - dxy * v[2] + dxy * v[3];
817 Av[1] = (dyy*Real(2.0) + beta[1]) * v[1] + dxy * v[2] - dxy * v[3];
818 Av[2] = -dxy * v[0] + dxy * v[1] + (dxx*Real(2.0) + beta[2]) * v[2];
819 Av[3] = dxy * v[0] - dxy * v[1] + (dxx*Real(2.0) + beta[3]) * v[3];
820 };
821 Real sol[4] = {0, 0, 0, 0};
822 pcg_solve<4>(sol, b.data(), mat, precond, 8, Real(1.e-8));
823 ex(i-1,j ,k ) = sol[0];
824 ex(i ,j ,k ) = sol[1];
825 ey(i ,j-1,k ) = sol[2];
826 ey(i ,j ,k ) = sol[3];
827 } else {
828 LUSolver<4,Real> lusolver
829 ({dyy*Real(2.0) + beta[0],
830 Real(0.0),
831 -dxy,
832 dxy,
833 //
834 Real(0.0),
835 dyy*Real(2.0) + beta[1],
836 dxy,
837 -dxy,
838 //
839 -dxy,
840 dxy,
841 dxx*Real(2.0) + beta[2],
842 Real(0.0),
843 //
844 dxy,
845 -dxy,
846 Real(0.0),
847 dxx*Real(2.0) + beta[3]});
848 lusolver(beta.data(), b.data());
849 ex(i-1,j ,k ) = beta[0];
850 ex(i ,j ,k ) = beta[1];
851 ey(i ,j-1,k ) = beta[2];
852 ey(i ,j ,k ) = beta[3];
853 }
854
855#elif (AMREX_SPACEDIM == 3)
856
857 if (my_color != color) { return; }
858
859 Real dxx = adxinv[0]*adxinv[0];
860 Real dyy = adxinv[1]*adxinv[1];
861 Real dzz = adxinv[2]*adxinv[2];
862 Real dxy = adxinv[0]*adxinv[1];
863 Real dxz = adxinv[0]*adxinv[2];
864 Real dyz = adxinv[1]*adxinv[2];
865
866 GpuArray<Real,6> b
867 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
868 ex(i-1,j+1,k ))
869 - dzz * ( ex(i-1,j ,k+1) +
870 ex(i-1,j ,k-1))
871 + dxy * ( ey(i-1,j-1,k )
872 -ey(i-1,j ,k ))
873 + dxz * ( ez(i-1,j ,k-1)
874 -ez(i-1,j ,k ))),
875 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
876 ex(i ,j+1,k ))
877 - dzz * ( ex(i ,j ,k+1) +
878 ex(i ,j ,k-1))
879 + dxy * (-ey(i+1,j-1,k )
880 +ey(i+1,j ,k ))
881 + dxz * (-ez(i+1,j ,k-1)
882 +ez(i+1,j ,k ))),
883 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
884 ey(i+1,j-1,k ))
885 - dzz * ( ey(i ,j-1,k-1) +
886 ey(i ,j-1,k+1))
887 + dxy * ( ex(i-1,j-1,k )
888 -ex(i ,j-1,k ))
889 + dyz * ( ez(i ,j-1,k-1)
890 -ez(i ,j-1,k ))),
891 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
892 ey(i+1,j ,k ))
893 - dzz * ( ey(i ,j ,k-1) +
894 ey(i ,j ,k+1))
895 + dxy * (-ex(i-1,j+1,k )
896 +ex(i ,j+1,k ))
897 + dyz * (-ez(i ,j+1,k-1)
898 +ez(i ,j+1,k ))),
899 rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) +
900 ez(i+1,j ,k-1))
901 - dyy * ( ez(i ,j-1,k-1) +
902 ez(i ,j+1,k-1))
903 + dxz * ( ex(i-1,j ,k-1)
904 -ex(i ,j ,k-1))
905 + dyz * ( ey(i ,j-1,k-1)
906 -ey(i ,j ,k-1))),
907 rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) +
908 ez(i+1,j ,k ))
909 - dyy * ( ez(i ,j-1,k ) +
910 ez(i ,j+1,k ))
911 + dxz * (-ex(i-1,j ,k+1)
912 +ex(i ,j ,k+1))
913 + dyz * (-ey(i ,j-1,k+1)
914 +ey(i ,j ,k+1)))};
915
916 GpuArray<Real,6> beta;
917
918 if (sinfo.xlo_is_symmetric(i)) {
919 b[0] = -b[1];
920 beta[0] = beta[1] = betax(i,j,k);
921 } else if (sinfo.xhi_is_symmetric(i)) {
922 b[1] = -b[0];
923 beta[0] = beta[1] = betax(i-1,j,k);
924 } else {
925 beta[0] = betax(i-1,j,k);
926 beta[1] = betax(i ,j,k);
927 }
928
929 if (sinfo.ylo_is_symmetric(j)) {
930 b[2] = -b[3];
931 beta[2] = beta[3] = betay(i,j,k);
932 } else if (sinfo.yhi_is_symmetric(j)) {
933 b[3] = -b[2];
934 beta[2] = beta[3] = betay(i,j-1,k);
935 } else {
936 beta[2] = betay(i,j-1,k);
937 beta[3] = betay(i,j ,k);
938 }
939
940 if (sinfo.zlo_is_symmetric(k)) {
941 b[4] = -b[5];
942 beta[4] = beta[5] = betaz(i,j,k);
943 } else if (sinfo.zhi_is_symmetric(k)) {
944 b[5] = -b[4];
945 beta[4] = beta[5] = betaz(i,j,k-1);
946 } else {
947 beta[4] = betaz(i,j,k-1);
948 beta[5] = betaz(i,j,k );
949 }
950
951 if constexpr (PCG) {
952 Real diagInv[6] = {Real(1.0) / ((dyy+dzz)*Real(2.0) + beta[0]),
953 Real(1.0) / ((dyy+dzz)*Real(2.0) + beta[1]),
954 Real(1.0) / ((dxx+dzz)*Real(2.0) + beta[2]),
955 Real(1.0) / ((dxx+dzz)*Real(2.0) + beta[3]),
956 Real(1.0) / ((dxx+dyy)*Real(2.0) + beta[4]),
957 Real(1.0) / ((dxx+dyy)*Real(2.0) + beta[5])};
958 auto precond = [&] (Real * AMREX_RESTRICT z,
959 Real const* AMREX_RESTRICT r)
960 {
961 for (int m = 0; m < 6; ++m) { z[m] = r[m] * diagInv[m]; }
962 };
963 auto mat = [&] (Real * AMREX_RESTRICT Av,
964 Real const* AMREX_RESTRICT v)
965 {
966 Av[0] = ((dyy+dzz)*Real(2.0) + beta[0]) * v[0] - dxy * v[2]
967 + dxy * v[3] - dxz * v[4] + dxz * v[5];
968 Av[1] = ((dyy+dzz)*Real(2.0) + beta[1]) * v[1] + dxy * v[2]
969 - dxy * v[3] + dxz * v[4] - dxz * v[5];
970 Av[2] = -dxy * v[0] + dxy * v[1] + ((dxx+dzz)*Real(2.0) + beta[2]) * v[2]
971 - dyz * v[4] + dyz * v[5];
972 Av[3] = dxy * v[0] - dxy * v[1] + ((dxx+dzz)*Real(2.0) + beta[3]) * v[3]
973 + dyz * v[4] - dyz * v[5];
974 Av[4] = -dxz * v[0] + dxz * v[1] - dyz * v[2] + dyz * v[3]
975 + ((dxx+dyy)*Real(2.0) + beta[4]) * v[4];
976 Av[5] = dxz * v[0] - dxz * v[1] + dyz * v[2] - dyz * v[3]
977 + ((dxx+dyy)*Real(2.0) + beta[5]) * v[5];
978 };
979 Real sol[6] = {0, 0, 0, 0, 0, 0};
980 pcg_solve<6>(sol, b.data(), mat, precond, 8, Real(1.e-8));
981 ex(i-1,j ,k ) = sol[0];
982 ex(i ,j ,k ) = sol[1];
983 ey(i ,j-1,k ) = sol[2];
984 ey(i ,j ,k ) = sol[3];
985 ez(i ,j ,k-1) = sol[4];
986 ez(i ,j ,k ) = sol[5];
987 } else {
988 LUSolver<6,Real> lusolver
989 ({(dyy+dzz)*Real(2.0) + beta[0],
990 Real(0.0),
991 -dxy,
992 dxy,
993 -dxz,
994 dxz,
995 //
996 Real(0.0),
997 (dyy+dzz)*Real(2.0) + beta[1],
998 dxy,
999 -dxy,
1000 dxz,
1001 -dxz,
1002 //
1003 -dxy,
1004 dxy,
1005 (dxx+dzz)*Real(2.0) + beta[2],
1006 Real(0.0),
1007 -dyz,
1008 dyz,
1009 //
1010 dxy,
1011 -dxy,
1012 Real(0.0),
1013 (dxx+dzz)*Real(2.0) + beta[3],
1014 dyz,
1015 -dyz,
1016 //
1017 -dxz,
1018 dxz,
1019 -dyz,
1020 dyz,
1021 (dxx+dyy)*Real(2.0) + beta[4],
1022 Real(0.0),
1023 //
1024 dxz,
1025 -dxz,
1026 dyz,
1027 -dyz,
1028 Real(0.0),
1029 (dxx+dyy)*Real(2.0) + beta[5]});
1030 lusolver(beta.data(), b.data());
1031 ex(i-1,j ,k ) = beta[0];
1032 ex(i ,j ,k ) = beta[1];
1033 ey(i ,j-1,k ) = beta[2];
1034 ey(i ,j ,k ) = beta[3];
1035 ez(i ,j ,k-1) = beta[4];
1036 ez(i ,j ,k ) = beta[5];
1037 }
1038#endif
1039}
1040
1041#endif
1042
1043AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1044void mlcurlcurl_interpadd (int dir, int i, int j, int k,
1045 Array4<Real> const& fine,
1046 Array4<Real const> const& crse)
1047{
1048 int ic = amrex::coarsen(i,2);
1049 int jc = amrex::coarsen(j,2);
1050 int kc = amrex::coarsen(k,2);
1051 if (dir == 0) {
1052 bool j_is_odd = (jc*2 != j);
1053 bool k_is_odd = (kc*2 != k);
1054 if (j_is_odd && k_is_odd) {
1055 fine(i,j,k) += Real(0.25) *
1056 (crse(ic,jc,kc ) + crse(ic,jc+1,kc ) +
1057 crse(ic,jc,kc+1) + crse(ic,jc+1,kc+1));
1058 } else if (j_is_odd) {
1059 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc));
1060 } else if (k_is_odd) {
1061 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1));
1062 } else {
1063 fine(i,j,k) += crse(ic,jc,kc);
1064 }
1065 } else if (dir == 1) {
1066 bool i_is_odd = (ic*2 != i);
1067 bool k_is_odd = (kc*2 != k);
1068 if (i_is_odd && k_is_odd) {
1069 fine(i,j,k) += Real(0.25) *
1070 (crse(ic ,jc,kc ) + crse(ic+1,jc,kc ) +
1071 crse(ic ,jc,kc+1) + crse(ic+1,jc,kc+1));
1072 } else if (i_is_odd) {
1073 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc));
1074 } else if (k_is_odd) {
1075 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1));
1076 } else {
1077 fine(i,j,k) += crse(ic,jc,kc);
1078 }
1079 } else {
1080 bool i_is_odd = (ic*2 != i);
1081 bool j_is_odd = (jc*2 != j);
1082 if (i_is_odd && j_is_odd) {
1083 fine(i,j,k) += Real(0.25) *
1084 (crse(ic ,jc ,kc) + crse(ic+1,jc ,kc) +
1085 crse(ic ,jc+1,kc) + crse(ic+1,jc+1,kc));
1086 } else if (i_is_odd) {
1087 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc));
1088 } else if (j_is_odd) {
1089 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc));
1090 } else {
1091 fine(i,j,k) += crse(ic,jc,kc);
1092 }
1093 }
1094}
1095
1096AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1097void mlcurlcurl_restriction (int dir, int i, int j, int k,
1098 Array4<Real> const& crse,
1099 Array4<Real const> const& fine,
1100 CurlCurlDirichletInfo const& dinfo)
1101{
1102 int ii = i*2;
1103 int jj = j*2;
1104 int kk = k*2;
1105 if (dinfo.is_dirichlet_edge(dir,ii,jj,kk)) {
1106 crse(i,j,k) = Real(0.0);
1107 }
1108 else
1109 {
1110#if (AMREX_SPACEDIM == 3)
1111 if (dir == 0) {
1112 crse(i,j,k) = Real(1./32.) * (fine(ii ,jj-1,kk-1) +
1113 fine(ii ,jj ,kk-1) * Real(2.0) +
1114 fine(ii ,jj+1,kk-1) +
1115 fine(ii ,jj-1,kk ) * Real(2.0) +
1116 fine(ii ,jj ,kk ) * Real(4.0) +
1117 fine(ii ,jj+1,kk ) * Real(2.0) +
1118 fine(ii ,jj-1,kk+1) +
1119 fine(ii ,jj ,kk+1) * Real(2.0) +
1120 fine(ii ,jj+1,kk+1) +
1121 fine(ii+1,jj-1,kk-1) +
1122 fine(ii+1,jj ,kk-1) * Real(2.0) +
1123 fine(ii+1,jj+1,kk-1) +
1124 fine(ii+1,jj-1,kk ) * Real(2.0) +
1125 fine(ii+1,jj ,kk ) * Real(4.0) +
1126 fine(ii+1,jj+1,kk ) * Real(2.0) +
1127 fine(ii+1,jj-1,kk+1) +
1128 fine(ii+1,jj ,kk+1) * Real(2.0) +
1129 fine(ii+1,jj+1,kk+1) );
1130 } else if (dir == 1) {
1131 crse(i,j,k) = Real(1./32.) * (fine(ii-1,jj ,kk-1) +
1132 fine(ii ,jj ,kk-1) * Real(2.0) +
1133 fine(ii+1,jj ,kk-1) +
1134 fine(ii-1,jj ,kk ) * Real(2.0) +
1135 fine(ii ,jj ,kk ) * Real(4.0) +
1136 fine(ii+1,jj ,kk ) * Real(2.0) +
1137 fine(ii-1,jj ,kk+1) +
1138 fine(ii ,jj ,kk+1) * Real(2.0) +
1139 fine(ii+1,jj ,kk+1) +
1140 fine(ii-1,jj+1,kk-1) +
1141 fine(ii ,jj+1,kk-1) * Real(2.0) +
1142 fine(ii+1,jj+1,kk-1) +
1143 fine(ii-1,jj+1,kk ) * Real(2.0) +
1144 fine(ii ,jj+1,kk ) * Real(4.0) +
1145 fine(ii+1,jj+1,kk ) * Real(2.0) +
1146 fine(ii-1,jj+1,kk+1) +
1147 fine(ii ,jj+1,kk+1) * Real(2.0) +
1148 fine(ii+1,jj+1,kk+1) );
1149 } else {
1150 crse(i,j,k) = Real(1./32.) * (fine(ii-1,jj-1,kk ) +
1151 fine(ii ,jj-1,kk ) * Real(2.0) +
1152 fine(ii+1,jj-1,kk ) +
1153 fine(ii-1,jj ,kk ) * Real(2.0) +
1154 fine(ii ,jj ,kk ) * Real(4.0) +
1155 fine(ii+1,jj ,kk ) * Real(2.0) +
1156 fine(ii-1,jj+1,kk ) +
1157 fine(ii ,jj+1,kk ) * Real(2.0) +
1158 fine(ii+1,jj+1,kk ) +
1159 fine(ii-1,jj-1,kk+1) +
1160 fine(ii ,jj-1,kk+1) * Real(2.0) +
1161 fine(ii+1,jj-1,kk+1) +
1162 fine(ii-1,jj ,kk+1) * Real(2.0) +
1163 fine(ii ,jj ,kk+1) * Real(4.0) +
1164 fine(ii+1,jj ,kk+1) * Real(2.0) +
1165 fine(ii-1,jj+1,kk+1) +
1166 fine(ii ,jj+1,kk+1) * Real(2.0) +
1167 fine(ii+1,jj+1,kk+1) );
1168 }
1169#elif (AMREX_SPACEDIM == 2)
1170 if (dir == 0) {
1171 crse(i,j,0) = Real(0.125) * (fine(ii ,jj-1,0) +
1172 fine(ii ,jj ,0) * Real(2.0) +
1173 fine(ii ,jj+1,0) +
1174 fine(ii+1,jj-1,0) +
1175 fine(ii+1,jj ,0) * Real(2.0) +
1176 fine(ii+1,jj+1,0) );
1177 } else if (dir == 1) {
1178 crse(i,j,0) = Real(0.125) * (fine(ii-1,jj ,0) +
1179 fine(ii ,jj ,0) * Real(2.0) +
1180 fine(ii+1,jj ,0) +
1181 fine(ii-1,jj+1,0) +
1182 fine(ii ,jj+1,0) * Real(2.0) +
1183 fine(ii+1,jj+1,0) );
1184 } else {
1185 crse(i,j,0) = Real(1./16.)*(fine(ii-1,jj-1,0) + Real(2.)*fine(ii ,jj-1,0) + fine(ii+1,jj-1,0)
1186 + Real(2.)*fine(ii-1,jj ,0) + Real(4.)*fine(ii ,jj ,0) + Real(2.)*fine(ii+1,jj ,0)
1187 + fine(ii-1,jj+1,0) + Real(2.)*fine(ii ,jj+1,0) + fine(ii+1,jj+1,0));
1188
1189 }
1190
1191#elif (AMREX_SPACEDIM == 1)
1192 if (dir == 0) {
1193 crse(i,0,0) = Real(0.5) * (fine(ii,0,0) + fine(ii+1,0,0));
1194 } else {
1195 crse(i,0,0) = Real(0.25) * (fine(ii-1,0,0) +
1196 fine(ii ,0,0) * Real(2.0) +
1197 fine(ii+1,0,0));
1198 }
1199#endif
1200 }
1201}
1202
1203AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1204void mlcurlcurl_bc_symmetry (int i, int j, int k, Orientation face, IndexType it,
1205 Array4<Real> const& a)
1206{
1207 int const idir = face.coordDir();
1208 int offset = face.isLow() ? 1 : -1;
1209 Real sign;
1210 if (it.cellCentered(idir)) {
1211 sign = Real(-1.0);
1212 } else {
1213 sign = Real(1.0);
1214 offset *= 2;
1215 }
1216
1217 if (idir == 0) {
1218 a(i,j,k) = sign * a(i+offset,j,k);
1219 } else if (idir == 1) {
1220 a(i,j,k) = sign * a(i,j+offset,k);
1221 } else {
1222 a(i,j,k) = sign * a(i,j,k+offset);
1223 }
1224}
1225
1226}
1227
1228#endif
AMReX_Array4.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMREX_RESTRICT
#define AMREX_RESTRICT
Definition AMReX_Extension.H:37
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_LUSolver.H
AMReX_PCGSolver.H
amrex::IndexTypeND::cellCentered
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool cellCentered() const noexcept
True if the IndexTypeND is CELL based in all directions.
Definition AMReX_IndexType.H:101
amrex::Orientation
Encapsulation of the Orientation of the Faces of a Box.
Definition AMReX_Orientation.H:29
amrex::Orientation::coordDir
AMREX_GPU_HOST_DEVICE int coordDir() const noexcept
Returns the coordinate direction.
Definition AMReX_Orientation.H:83
amrex::Orientation::isLow
AMREX_GPU_HOST_DEVICE bool isLow() const noexcept
Returns true if Orientation is low.
Definition AMReX_Orientation.H:89
amrex
Definition AMReX_Amr.cpp:49
amrex::CurlCurlStateType
CurlCurlStateType
Definition AMReX_MLCurlCurl_K.H:326
amrex::mlcurlcurl_interpadd
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_interpadd(int dir, int i, int j, int k, Array4< Real > const &fine, Array4< Real const > const &crse)
Definition AMReX_MLCurlCurl_K.H:1044
amrex::mlcurlcurl_1D
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_1D(int i, int j, int k, Array4< Real > const &ex, Array4< Real > const &ey, Array4< Real > const &ez, Array4< Real const > const &rhsx, Array4< Real const > const &rhsy, Array4< Real const > const &rhsz, Real beta, GpuArray< Real, AMREX_SPACEDIM > const &adxinv, int color, CurlCurlDirichletInfo const &dinfo)
Definition AMReX_MLCurlCurl_K.H:463
amrex::coarsen
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo/rati...
Definition AMReX_Box.H:1304
amrex::mlcurlcurl_adotx_y
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_adotx_y(int i, int j, int k, Array4< Real > const &Ay, Array4< Real const > const &ex, Array4< Real const > const &ey, Array4< Real const > const &ez, Real beta, GpuArray< Real, AMREX_SPACEDIM > const &adxinv)
Definition AMReX_MLCurlCurl_K.H:372
amrex::ignore_unused
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:111
amrex::mlcurlcurl_bc_symmetry
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_bc_symmetry(int i, int j, int k, Orientation face, IndexType it, Array4< Real > const &a)
Definition AMReX_MLCurlCurl_K.H:1204
amrex::mlcurlcurl_adotx_z
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_adotx_z(int i, int j, int k, Array4< Real > const &Az, Array4< Real const > const &ex, Array4< Real const > const &ey, Array4< Real const > const &ez, Real beta, GpuArray< Real, AMREX_SPACEDIM > const &adxinv)
Definition AMReX_MLCurlCurl_K.H:418
amrex::mlcurlcurl_adotx_x
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_adotx_x(int i, int j, int k, Array4< Real > const &Ax, Array4< Real const > const &ex, Array4< Real const > const &ey, Array4< Real const > const &ez, Real beta, GpuArray< Real, AMREX_SPACEDIM > const &adxinv)
Definition AMReX_MLCurlCurl_K.H:329
amrex::mlcurlcurl_restriction
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_restriction(int dir, int i, int j, int k, Array4< Real > const &crse, Array4< Real const > const &fine, CurlCurlDirichletInfo const &dinfo)
Definition AMReX_MLCurlCurl_K.H:1097
amrex::Array4
Definition AMReX_Array4.H:61
amrex::CurlCurlDirichletInfo
Definition AMReX_MLCurlCurl_K.H:185
amrex::CurlCurlDirichletInfo::is_dirichlet_node
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_node(int i, int j, int k) const
Definition AMReX_MLCurlCurl_K.H:190
amrex::CurlCurlDirichletInfo::is_dirichlet_edge
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_edge(int dim, int i, int j, int k) const
Definition AMReX_MLCurlCurl_K.H:246
amrex::CurlCurlDirichletInfo::is_dirichlet_z_edge
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_z_edge(int i, int j, int) const
Definition AMReX_MLCurlCurl_K.H:234
amrex::CurlCurlDirichletInfo::dirichlet_hi
IntVect dirichlet_hi
Definition AMReX_MLCurlCurl_K.H:187
amrex::CurlCurlDirichletInfo::is_dirichlet_y_edge
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_y_edge(int i, int, int k) const
Definition AMReX_MLCurlCurl_K.H:222
amrex::CurlCurlDirichletInfo::dirichlet_lo
IntVect dirichlet_lo
Definition AMReX_MLCurlCurl_K.H:186
amrex::CurlCurlDirichletInfo::is_dirichlet_x_edge
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_x_edge(int, int j, int k) const
Definition AMReX_MLCurlCurl_K.H:207
amrex::CurlCurlSymmetryInfo
Definition AMReX_MLCurlCurl_K.H:259
amrex::CurlCurlSymmetryInfo::xlo_is_symmetric
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool xlo_is_symmetric(int i) const
Definition AMReX_MLCurlCurl_K.H:264
amrex::CurlCurlSymmetryInfo::xhi_is_symmetric
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool xhi_is_symmetric(int i) const
Definition AMReX_MLCurlCurl_K.H:270
amrex::CurlCurlSymmetryInfo::symmetry_hi
IntVect symmetry_hi
Definition AMReX_MLCurlCurl_K.H:261
amrex::CurlCurlSymmetryInfo::symmetry_lo
IntVect symmetry_lo
Definition AMReX_MLCurlCurl_K.H:260
amrex::CurlCurlSymmetryInfo::is_symmetric
bool is_symmetric(int dir, int side, int idx) const
Definition AMReX_MLCurlCurl_K.H:303
amrex::GpuArray
Definition AMReX_Array.H:34