Block-Structured AMR Software Framework
 
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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, bool valid_x)
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 if (valid_x) {
484 ex(i,j,k) = rhsx(i,j,k) / beta;
485 }
486
487 Real gamma_y = dxx * Real(2.0) + beta;
488 Real ccey = - dxx * (ey(i-1,j ,k ) +
489 ey(i+1,j ,k ) );
490 Real res_y = rhsy(i,j,k) - ( gamma_y * ey(i,j,k) + ccey );
491 ey(i,j,k) += res_y/gamma_y;
492
493 Real gamma_z = dxx * Real(2.0) + beta;
494 Real ccez = -dxx * (ez(i-1,j ,k ) +
495 ez(i+1,j ,k ) );
496 Real res_z = rhsz(i,j,k) - ( gamma_z * ez(i,j,k) + ccez );
497 ez(i,j,k) += res_z/gamma_z;
498 }
499}
500
501
502AMREX_GPU_DEVICE AMREX_FORCE_INLINE
503void mlcurlcurl_1D (int i, int j, int k,
504 Array4<Real> const& ex,
505 Array4<Real> const& ey,
506 Array4<Real> const& ez,
507 Array4<Real const> const& rhsx,
508 Array4<Real const> const& rhsy,
509 Array4<Real const> const& rhsz,
510 Array4<Real const> const& betax,
511 Array4<Real const> const& betay,
512 Array4<Real const> const& betaz,
513 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
514 int color,
515 CurlCurlDirichletInfo const& dinfo, bool valid_x)
516{
517 if (dinfo.is_dirichlet_node(i,j,k)) {return; }
518 Real dxx = adxinv[0] * adxinv[0];
519
520 int my_color = i%2;
521
522 if (my_color == color)
523 {
524 if (valid_x) {
525 ex(i,j,k) = rhsx(i,j,k) / betax(i,j,k);
526 }
527
528 Real gamma_y = dxx * Real(2.0) + betay(i,j,k);
529 Real ccey = - dxx * (ey(i-1,j ,k ) +
530 ey(i+1,j ,k ) );
531 Real res_y = rhsy(i,j,k) - ( gamma_y * ey(i,j,k) + ccey );
532 ey(i,j,k) += res_y/gamma_y;
533
534 Real gamma_z = dxx * Real(2.0) + betaz(i,j,k);
535 Real ccez = -dxx * (ez(i-1,j ,k ) +
536 ez(i+1,j ,k ) );
537 Real res_z = rhsz(i,j,k) - ( gamma_z * ez(i,j,k) + ccez );
538 ez(i,j,k) += res_z/gamma_z;
539 }
540}
541
542
543#endif
544
545#if (AMREX_SPACEDIM > 1)
546
547AMREX_GPU_DEVICE AMREX_FORCE_INLINE
548void mlcurlcurl_gs4_lu (int i, int j, int k,
549 Array4<Real> const& ex,
550 Array4<Real> const& ey,
551 Array4<Real> const& ez,
552 Array4<Real const> const& rhsx,
553 Array4<Real const> const& rhsy,
554 Array4<Real const> const& rhsz,
555#if (AMREX_SPACEDIM == 2)
556 Real beta,
557#endif
558 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
559 int color, LUSolver<AMREX_SPACEDIM*2,Real> const& lusolver,
560 CurlCurlDirichletInfo const& dinfo,
561 CurlCurlSymmetryInfo const& sinfo)
562{
563 if (dinfo.is_dirichlet_node(i,j,k)) { return; }
564
565 int my_color = i%2 + 2*(j%2);
566 if (k%2 != 0) {
567 my_color = 3 - my_color;
568 }
569
570#if (AMREX_SPACEDIM == 2)
571
572 Real dxx = adxinv[0] * adxinv[0];
573 Real dyy = adxinv[1] * adxinv[1];
574
575 if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) ||
576 ((my_color == 1 || my_color == 2) && (color == 1 || color == 2)))
577 {
578 Real gamma = (dxx+dyy)*Real(2.0) + beta;
579 Real ccez = - dxx * (ez(i-1,j ,k ) +
580 ez(i+1,j ,k ))
581 - dyy * (ez(i ,j-1,k ) +
582 ez(i ,j+1,k ));
583 Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez);
584 constexpr Real omega = Real(1.15);
585 ez(i,j,k) += omega/gamma * res;
586 }
587
588 if (my_color != color) { return; }
589
590 Real dxy = adxinv[0]*adxinv[1];
591
592 GpuArray<Real,6> b
593 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
594 ex(i-1,j+1,k ))
595 + dxy * ( ey(i-1,j-1,k )
596 -ey(i-1,j ,k ))),
597 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
598 ex(i ,j+1,k ))
599 + dxy * (-ey(i+1,j-1,k )
600 +ey(i+1,j ,k ))),
601 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
602 ey(i+1,j-1,k ))
603 + dxy * ( ex(i-1,j-1,k )
604 -ex(i ,j-1,k ))),
605 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
606 ey(i+1,j ,k ))
607 + dxy * (-ex(i-1,j+1,k )
608 +ex(i ,j+1,k )))};
609
610 if (sinfo.xlo_is_symmetric(i)) {
611 b[0] = -b[1];
612 } else if (sinfo.xhi_is_symmetric(i)) {
613 b[1] = -b[0];
614 }
615
616 if (sinfo.ylo_is_symmetric(j)) {
617 b[2] = -b[3];
618 } else if (sinfo.yhi_is_symmetric(j)) {
619 b[3] = -b[2];
620 }
621
622 GpuArray<Real,4> x;
623 lusolver(x.data(), b.data());
624 ex(i-1,j ,k ) = x[0];
625 ex(i ,j ,k ) = x[1];
626 ey(i ,j-1,k ) = x[2];
627 ey(i ,j ,k ) = x[3];
628
629#elif (AMREX_SPACEDIM == 3)
630
631 if (my_color != color) { return; }
632
633 Real dxx = adxinv[0]*adxinv[0];
634 Real dyy = adxinv[1]*adxinv[1];
635 Real dzz = adxinv[2]*adxinv[2];
636 Real dxy = adxinv[0]*adxinv[1];
637 Real dxz = adxinv[0]*adxinv[2];
638 Real dyz = adxinv[1]*adxinv[2];
639
640 GpuArray<Real,6> b
641 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
642 ex(i-1,j+1,k ))
643 - dzz * ( ex(i-1,j ,k+1) +
644 ex(i-1,j ,k-1))
645 + dxy * ( ey(i-1,j-1,k )
646 -ey(i-1,j ,k ))
647 + dxz * ( ez(i-1,j ,k-1)
648 -ez(i-1,j ,k ))),
649 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
650 ex(i ,j+1,k ))
651 - dzz * ( ex(i ,j ,k+1) +
652 ex(i ,j ,k-1))
653 + dxy * (-ey(i+1,j-1,k )
654 +ey(i+1,j ,k ))
655 + dxz * (-ez(i+1,j ,k-1)
656 +ez(i+1,j ,k ))),
657 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
658 ey(i+1,j-1,k ))
659 - dzz * ( ey(i ,j-1,k-1) +
660 ey(i ,j-1,k+1))
661 + dxy * ( ex(i-1,j-1,k )
662 -ex(i ,j-1,k ))
663 + dyz * ( ez(i ,j-1,k-1)
664 -ez(i ,j-1,k ))),
665 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
666 ey(i+1,j ,k ))
667 - dzz * ( ey(i ,j ,k-1) +
668 ey(i ,j ,k+1))
669 + dxy * (-ex(i-1,j+1,k )
670 +ex(i ,j+1,k ))
671 + dyz * (-ez(i ,j+1,k-1)
672 +ez(i ,j+1,k ))),
673 rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) +
674 ez(i+1,j ,k-1))
675 - dyy * ( ez(i ,j-1,k-1) +
676 ez(i ,j+1,k-1))
677 + dxz * ( ex(i-1,j ,k-1)
678 -ex(i ,j ,k-1))
679 + dyz * ( ey(i ,j-1,k-1)
680 -ey(i ,j ,k-1))),
681 rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) +
682 ez(i+1,j ,k ))
683 - dyy * ( ez(i ,j-1,k ) +
684 ez(i ,j+1,k ))
685 + dxz * (-ex(i-1,j ,k+1)
686 +ex(i ,j ,k+1))
687 + dyz * (-ey(i ,j-1,k+1)
688 +ey(i ,j ,k+1)))};
689
690 if (sinfo.xlo_is_symmetric(i)) {
691 b[0] = -b[1];
692 } else if (sinfo.xhi_is_symmetric(i)) {
693 b[1] = -b[0];
694 }
695
696 if (sinfo.ylo_is_symmetric(j)) {
697 b[2] = -b[3];
698 } else if (sinfo.yhi_is_symmetric(j)) {
699 b[3] = -b[2];
700 }
701
702 if (sinfo.zlo_is_symmetric(k)) {
703 b[4] = -b[5];
704 } else if (sinfo.zhi_is_symmetric(k)) {
705 b[5] = -b[4];
706 }
707
708 GpuArray<Real,6> x;
709 lusolver(x.data(), b.data());
710 ex(i-1,j ,k ) = x[0];
711 ex(i ,j ,k ) = x[1];
712 ey(i ,j-1,k ) = x[2];
713 ey(i ,j ,k ) = x[3];
714 ez(i ,j ,k-1) = x[4];
715 ez(i ,j ,k ) = x[5];
716#endif
717}
718
719template <bool PCG>
720AMREX_GPU_DEVICE AMREX_FORCE_INLINE
721void mlcurlcurl_gs4 (int i, int j, int k,
722 Array4<Real> const& ex,
723 Array4<Real> const& ey,
724 Array4<Real> const& ez,
725 Array4<Real const> const& rhsx,
726 Array4<Real const> const& rhsy,
727 Array4<Real const> const& rhsz,
728 GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
729 int color,
730 Array4<Real const> const& betax,
731 Array4<Real const> const& betay,
732 Array4<Real const> const& betaz,
733 CurlCurlDirichletInfo const& dinfo,
734 CurlCurlSymmetryInfo const& sinfo)
735{
736 if (dinfo.is_dirichlet_node(i,j,k)) { return; }
737
738 int my_color = i%2 + 2*(j%2);
739 if (k%2 != 0) {
740 my_color = 3 - my_color;
741 }
742
743#if (AMREX_SPACEDIM == 2)
744
745 Real dxx = adxinv[0] * adxinv[0];
746 Real dyy = adxinv[1] * adxinv[1];
747
748 if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) ||
749 ((my_color == 1 || my_color == 2) && (color == 1 || color == 2)))
750 {
751 Real gamma = (dxx+dyy)*Real(2.0) + betaz(i,j,k);
752 Real ccez = - dxx * (ez(i-1,j ,k ) +
753 ez(i+1,j ,k ))
754 - dyy * (ez(i ,j-1,k ) +
755 ez(i ,j+1,k ));
756 Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez);
757 constexpr Real omega = Real(1.15);
758 ez(i,j,k) += omega/gamma * res;
759 }
760
761 if (my_color != color) { return; }
762
763 Real dxy = adxinv[0]*adxinv[1];
764
765 GpuArray<Real,6> b
766 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
767 ex(i-1,j+1,k ))
768 + dxy * ( ey(i-1,j-1,k )
769 -ey(i-1,j ,k ))),
770 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
771 ex(i ,j+1,k ))
772 + dxy * (-ey(i+1,j-1,k )
773 +ey(i+1,j ,k ))),
774 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
775 ey(i+1,j-1,k ))
776 + dxy * ( ex(i-1,j-1,k )
777 -ex(i ,j-1,k ))),
778 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
779 ey(i+1,j ,k ))
780 + dxy * (-ex(i-1,j+1,k )
781 +ex(i ,j+1,k )))};
782
783 GpuArray<Real,4> beta;
784
785 if (sinfo.xlo_is_symmetric(i)) {
786 b[0] = -b[1];
787 beta[0] = beta[1] = betax(i,j,k);
788 } else if (sinfo.xhi_is_symmetric(i)) {
789 b[1] = -b[0];
790 beta[0] = beta[1] = betax(i-1,j,k);
791 } else {
792 beta[0] = betax(i-1,j,k);
793 beta[1] = betax(i ,j,k);
794 }
795
796 if (sinfo.ylo_is_symmetric(j)) {
797 b[2] = -b[3];
798 beta[2] = beta[3] = betay(i,j,k);
799 } else if (sinfo.yhi_is_symmetric(j)) {
800 b[3] = -b[2];
801 beta[2] = beta[3] = betay(i,j-1,k);
802 } else {
803 beta[2] = betay(i,j-1,k);
804 beta[3] = betay(i,j ,k);
805 }
806
807 if constexpr (PCG) {
808 Real diagInv[4] = {Real(1.0) / (dyy*Real(2.0) + beta[0]),
809 Real(1.0) / (dyy*Real(2.0) + beta[1]),
810 Real(1.0) / (dxx*Real(2.0) + beta[2]),
811 Real(1.0) / (dxx*Real(2.0) + beta[3])};
812 auto precond = [&] (Real * AMREX_RESTRICT z,
813 Real const* AMREX_RESTRICT r)
814 {
815 for (int m = 0; m < 4; ++m) { z[m] = r[m] * diagInv[m]; }
816 };
817 auto mat = [&] (Real * AMREX_RESTRICT Av,
818 Real const* AMREX_RESTRICT v)
819 {
820 Av[0] = (dyy*Real(2.0) + beta[0]) * v[0] - dxy * v[2] + dxy * v[3];
821 Av[1] = (dyy*Real(2.0) + beta[1]) * v[1] + dxy * v[2] - dxy * v[3];
822 Av[2] = -dxy * v[0] + dxy * v[1] + (dxx*Real(2.0) + beta[2]) * v[2];
823 Av[3] = dxy * v[0] - dxy * v[1] + (dxx*Real(2.0) + beta[3]) * v[3];
824 };
825 Real sol[4] = {0, 0, 0, 0};
826 pcg_solve<4>(sol, b.data(), mat, precond, 8, Real(1.e-8));
827 ex(i-1,j ,k ) = sol[0];
828 ex(i ,j ,k ) = sol[1];
829 ey(i ,j-1,k ) = sol[2];
830 ey(i ,j ,k ) = sol[3];
831 } else {
832 LUSolver<4,Real> lusolver
833 ({dyy*Real(2.0) + beta[0],
834 Real(0.0),
835 -dxy,
836 dxy,
837 //
838 Real(0.0),
839 dyy*Real(2.0) + beta[1],
840 dxy,
841 -dxy,
842 //
843 -dxy,
844 dxy,
845 dxx*Real(2.0) + beta[2],
846 Real(0.0),
847 //
848 dxy,
849 -dxy,
850 Real(0.0),
851 dxx*Real(2.0) + beta[3]});
852 lusolver(beta.data(), b.data());
853 ex(i-1,j ,k ) = beta[0];
854 ex(i ,j ,k ) = beta[1];
855 ey(i ,j-1,k ) = beta[2];
856 ey(i ,j ,k ) = beta[3];
857 }
858
859#elif (AMREX_SPACEDIM == 3)
860
861 if (my_color != color) { return; }
862
863 Real dxx = adxinv[0]*adxinv[0];
864 Real dyy = adxinv[1]*adxinv[1];
865 Real dzz = adxinv[2]*adxinv[2];
866 Real dxy = adxinv[0]*adxinv[1];
867 Real dxz = adxinv[0]*adxinv[2];
868 Real dyz = adxinv[1]*adxinv[2];
869
870 GpuArray<Real,6> b
871 {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
872 ex(i-1,j+1,k ))
873 - dzz * ( ex(i-1,j ,k+1) +
874 ex(i-1,j ,k-1))
875 + dxy * ( ey(i-1,j-1,k )
876 -ey(i-1,j ,k ))
877 + dxz * ( ez(i-1,j ,k-1)
878 -ez(i-1,j ,k ))),
879 rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
880 ex(i ,j+1,k ))
881 - dzz * ( ex(i ,j ,k+1) +
882 ex(i ,j ,k-1))
883 + dxy * (-ey(i+1,j-1,k )
884 +ey(i+1,j ,k ))
885 + dxz * (-ez(i+1,j ,k-1)
886 +ez(i+1,j ,k ))),
887 rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
888 ey(i+1,j-1,k ))
889 - dzz * ( ey(i ,j-1,k-1) +
890 ey(i ,j-1,k+1))
891 + dxy * ( ex(i-1,j-1,k )
892 -ex(i ,j-1,k ))
893 + dyz * ( ez(i ,j-1,k-1)
894 -ez(i ,j-1,k ))),
895 rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
896 ey(i+1,j ,k ))
897 - dzz * ( ey(i ,j ,k-1) +
898 ey(i ,j ,k+1))
899 + dxy * (-ex(i-1,j+1,k )
900 +ex(i ,j+1,k ))
901 + dyz * (-ez(i ,j+1,k-1)
902 +ez(i ,j+1,k ))),
903 rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) +
904 ez(i+1,j ,k-1))
905 - dyy * ( ez(i ,j-1,k-1) +
906 ez(i ,j+1,k-1))
907 + dxz * ( ex(i-1,j ,k-1)
908 -ex(i ,j ,k-1))
909 + dyz * ( ey(i ,j-1,k-1)
910 -ey(i ,j ,k-1))),
911 rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) +
912 ez(i+1,j ,k ))
913 - dyy * ( ez(i ,j-1,k ) +
914 ez(i ,j+1,k ))
915 + dxz * (-ex(i-1,j ,k+1)
916 +ex(i ,j ,k+1))
917 + dyz * (-ey(i ,j-1,k+1)
918 +ey(i ,j ,k+1)))};
919
920 GpuArray<Real,6> beta;
921
922 if (sinfo.xlo_is_symmetric(i)) {
923 b[0] = -b[1];
924 beta[0] = beta[1] = betax(i,j,k);
925 } else if (sinfo.xhi_is_symmetric(i)) {
926 b[1] = -b[0];
927 beta[0] = beta[1] = betax(i-1,j,k);
928 } else {
929 beta[0] = betax(i-1,j,k);
930 beta[1] = betax(i ,j,k);
931 }
932
933 if (sinfo.ylo_is_symmetric(j)) {
934 b[2] = -b[3];
935 beta[2] = beta[3] = betay(i,j,k);
936 } else if (sinfo.yhi_is_symmetric(j)) {
937 b[3] = -b[2];
938 beta[2] = beta[3] = betay(i,j-1,k);
939 } else {
940 beta[2] = betay(i,j-1,k);
941 beta[3] = betay(i,j ,k);
942 }
943
944 if (sinfo.zlo_is_symmetric(k)) {
945 b[4] = -b[5];
946 beta[4] = beta[5] = betaz(i,j,k);
947 } else if (sinfo.zhi_is_symmetric(k)) {
948 b[5] = -b[4];
949 beta[4] = beta[5] = betaz(i,j,k-1);
950 } else {
951 beta[4] = betaz(i,j,k-1);
952 beta[5] = betaz(i,j,k );
953 }
954
955 if constexpr (PCG) {
956 Real diagInv[6] = {Real(1.0) / ((dyy+dzz)*Real(2.0) + beta[0]),
957 Real(1.0) / ((dyy+dzz)*Real(2.0) + beta[1]),
958 Real(1.0) / ((dxx+dzz)*Real(2.0) + beta[2]),
959 Real(1.0) / ((dxx+dzz)*Real(2.0) + beta[3]),
960 Real(1.0) / ((dxx+dyy)*Real(2.0) + beta[4]),
961 Real(1.0) / ((dxx+dyy)*Real(2.0) + beta[5])};
962 auto precond = [&] (Real * AMREX_RESTRICT z,
963 Real const* AMREX_RESTRICT r)
964 {
965 for (int m = 0; m < 6; ++m) { z[m] = r[m] * diagInv[m]; }
966 };
967 auto mat = [&] (Real * AMREX_RESTRICT Av,
968 Real const* AMREX_RESTRICT v)
969 {
970 Av[0] = ((dyy+dzz)*Real(2.0) + beta[0]) * v[0] - dxy * v[2]
971 + dxy * v[3] - dxz * v[4] + dxz * v[5];
972 Av[1] = ((dyy+dzz)*Real(2.0) + beta[1]) * v[1] + dxy * v[2]
973 - dxy * v[3] + dxz * v[4] - dxz * v[5];
974 Av[2] = -dxy * v[0] + dxy * v[1] + ((dxx+dzz)*Real(2.0) + beta[2]) * v[2]
975 - dyz * v[4] + dyz * v[5];
976 Av[3] = dxy * v[0] - dxy * v[1] + ((dxx+dzz)*Real(2.0) + beta[3]) * v[3]
977 + dyz * v[4] - dyz * v[5];
978 Av[4] = -dxz * v[0] + dxz * v[1] - dyz * v[2] + dyz * v[3]
979 + ((dxx+dyy)*Real(2.0) + beta[4]) * v[4];
980 Av[5] = dxz * v[0] - dxz * v[1] + dyz * v[2] - dyz * v[3]
981 + ((dxx+dyy)*Real(2.0) + beta[5]) * v[5];
982 };
983 Real sol[6] = {0, 0, 0, 0, 0, 0};
984 pcg_solve<6>(sol, b.data(), mat, precond, 8, Real(1.e-8));
985 ex(i-1,j ,k ) = sol[0];
986 ex(i ,j ,k ) = sol[1];
987 ey(i ,j-1,k ) = sol[2];
988 ey(i ,j ,k ) = sol[3];
989 ez(i ,j ,k-1) = sol[4];
990 ez(i ,j ,k ) = sol[5];
991 } else {
992 LUSolver<6,Real> lusolver
993 ({(dyy+dzz)*Real(2.0) + beta[0],
994 Real(0.0),
995 -dxy,
996 dxy,
997 -dxz,
998 dxz,
999 //
1000 Real(0.0),
1001 (dyy+dzz)*Real(2.0) + beta[1],
1002 dxy,
1003 -dxy,
1004 dxz,
1005 -dxz,
1006 //
1007 -dxy,
1008 dxy,
1009 (dxx+dzz)*Real(2.0) + beta[2],
1010 Real(0.0),
1011 -dyz,
1012 dyz,
1013 //
1014 dxy,
1015 -dxy,
1016 Real(0.0),
1017 (dxx+dzz)*Real(2.0) + beta[3],
1018 dyz,
1019 -dyz,
1020 //
1021 -dxz,
1022 dxz,
1023 -dyz,
1024 dyz,
1025 (dxx+dyy)*Real(2.0) + beta[4],
1026 Real(0.0),
1027 //
1028 dxz,
1029 -dxz,
1030 dyz,
1031 -dyz,
1032 Real(0.0),
1033 (dxx+dyy)*Real(2.0) + beta[5]});
1034 lusolver(beta.data(), b.data());
1035 ex(i-1,j ,k ) = beta[0];
1036 ex(i ,j ,k ) = beta[1];
1037 ey(i ,j-1,k ) = beta[2];
1038 ey(i ,j ,k ) = beta[3];
1039 ez(i ,j ,k-1) = beta[4];
1040 ez(i ,j ,k ) = beta[5];
1041 }
1042#endif
1043}
1044
1045#endif
1046
1047AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1048void mlcurlcurl_interpadd (int dir, int i, int j, int k,
1049 Array4<Real> const& fine,
1050 Array4<Real const> const& crse)
1051{
1052 int ic = amrex::coarsen(i,2);
1053 int jc = amrex::coarsen(j,2);
1054 int kc = amrex::coarsen(k,2);
1055 if (dir == 0) {
1056 bool j_is_odd = (jc*2 != j);
1057 bool k_is_odd = (kc*2 != k);
1058 if (j_is_odd && k_is_odd) {
1059 fine(i,j,k) += Real(0.25) *
1060 (crse(ic,jc,kc ) + crse(ic,jc+1,kc ) +
1061 crse(ic,jc,kc+1) + crse(ic,jc+1,kc+1));
1062 } else if (j_is_odd) {
1063 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc));
1064 } else if (k_is_odd) {
1065 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1));
1066 } else {
1067 fine(i,j,k) += crse(ic,jc,kc);
1068 }
1069 } else if (dir == 1) {
1070 bool i_is_odd = (ic*2 != i);
1071 bool k_is_odd = (kc*2 != k);
1072 if (i_is_odd && k_is_odd) {
1073 fine(i,j,k) += Real(0.25) *
1074 (crse(ic ,jc,kc ) + crse(ic+1,jc,kc ) +
1075 crse(ic ,jc,kc+1) + crse(ic+1,jc,kc+1));
1076 } else if (i_is_odd) {
1077 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc));
1078 } else if (k_is_odd) {
1079 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1));
1080 } else {
1081 fine(i,j,k) += crse(ic,jc,kc);
1082 }
1083 } else {
1084 bool i_is_odd = (ic*2 != i);
1085 bool j_is_odd = (jc*2 != j);
1086 if (i_is_odd && j_is_odd) {
1087 fine(i,j,k) += Real(0.25) *
1088 (crse(ic ,jc ,kc) + crse(ic+1,jc ,kc) +
1089 crse(ic ,jc+1,kc) + crse(ic+1,jc+1,kc));
1090 } else if (i_is_odd) {
1091 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc));
1092 } else if (j_is_odd) {
1093 fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc));
1094 } else {
1095 fine(i,j,k) += crse(ic,jc,kc);
1096 }
1097 }
1098}
1099
1100AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1101void mlcurlcurl_restriction (int dir, int i, int j, int k,
1102 Array4<Real> const& crse,
1103 Array4<Real const> const& fine,
1104 CurlCurlDirichletInfo const& dinfo)
1105{
1106 int ii = i*2;
1107 int jj = j*2;
1108 int kk = k*2;
1109 if (dinfo.is_dirichlet_edge(dir,ii,jj,kk)) {
1110 crse(i,j,k) = Real(0.0);
1111 }
1112 else
1113 {
1114#if (AMREX_SPACEDIM == 3)
1115 if (dir == 0) {
1116 crse(i,j,k) = Real(1./32.) * (fine(ii ,jj-1,kk-1) +
1117 fine(ii ,jj ,kk-1) * Real(2.0) +
1118 fine(ii ,jj+1,kk-1) +
1119 fine(ii ,jj-1,kk ) * Real(2.0) +
1120 fine(ii ,jj ,kk ) * Real(4.0) +
1121 fine(ii ,jj+1,kk ) * Real(2.0) +
1122 fine(ii ,jj-1,kk+1) +
1123 fine(ii ,jj ,kk+1) * Real(2.0) +
1124 fine(ii ,jj+1,kk+1) +
1125 fine(ii+1,jj-1,kk-1) +
1126 fine(ii+1,jj ,kk-1) * Real(2.0) +
1127 fine(ii+1,jj+1,kk-1) +
1128 fine(ii+1,jj-1,kk ) * Real(2.0) +
1129 fine(ii+1,jj ,kk ) * Real(4.0) +
1130 fine(ii+1,jj+1,kk ) * Real(2.0) +
1131 fine(ii+1,jj-1,kk+1) +
1132 fine(ii+1,jj ,kk+1) * Real(2.0) +
1133 fine(ii+1,jj+1,kk+1) );
1134 } else if (dir == 1) {
1135 crse(i,j,k) = Real(1./32.) * (fine(ii-1,jj ,kk-1) +
1136 fine(ii ,jj ,kk-1) * Real(2.0) +
1137 fine(ii+1,jj ,kk-1) +
1138 fine(ii-1,jj ,kk ) * Real(2.0) +
1139 fine(ii ,jj ,kk ) * Real(4.0) +
1140 fine(ii+1,jj ,kk ) * Real(2.0) +
1141 fine(ii-1,jj ,kk+1) +
1142 fine(ii ,jj ,kk+1) * Real(2.0) +
1143 fine(ii+1,jj ,kk+1) +
1144 fine(ii-1,jj+1,kk-1) +
1145 fine(ii ,jj+1,kk-1) * Real(2.0) +
1146 fine(ii+1,jj+1,kk-1) +
1147 fine(ii-1,jj+1,kk ) * Real(2.0) +
1148 fine(ii ,jj+1,kk ) * Real(4.0) +
1149 fine(ii+1,jj+1,kk ) * Real(2.0) +
1150 fine(ii-1,jj+1,kk+1) +
1151 fine(ii ,jj+1,kk+1) * Real(2.0) +
1152 fine(ii+1,jj+1,kk+1) );
1153 } else {
1154 crse(i,j,k) = Real(1./32.) * (fine(ii-1,jj-1,kk ) +
1155 fine(ii ,jj-1,kk ) * Real(2.0) +
1156 fine(ii+1,jj-1,kk ) +
1157 fine(ii-1,jj ,kk ) * Real(2.0) +
1158 fine(ii ,jj ,kk ) * Real(4.0) +
1159 fine(ii+1,jj ,kk ) * Real(2.0) +
1160 fine(ii-1,jj+1,kk ) +
1161 fine(ii ,jj+1,kk ) * Real(2.0) +
1162 fine(ii+1,jj+1,kk ) +
1163 fine(ii-1,jj-1,kk+1) +
1164 fine(ii ,jj-1,kk+1) * Real(2.0) +
1165 fine(ii+1,jj-1,kk+1) +
1166 fine(ii-1,jj ,kk+1) * Real(2.0) +
1167 fine(ii ,jj ,kk+1) * Real(4.0) +
1168 fine(ii+1,jj ,kk+1) * Real(2.0) +
1169 fine(ii-1,jj+1,kk+1) +
1170 fine(ii ,jj+1,kk+1) * Real(2.0) +
1171 fine(ii+1,jj+1,kk+1) );
1172 }
1173#elif (AMREX_SPACEDIM == 2)
1174 if (dir == 0) {
1175 crse(i,j,0) = Real(0.125) * (fine(ii ,jj-1,0) +
1176 fine(ii ,jj ,0) * Real(2.0) +
1177 fine(ii ,jj+1,0) +
1178 fine(ii+1,jj-1,0) +
1179 fine(ii+1,jj ,0) * Real(2.0) +
1180 fine(ii+1,jj+1,0) );
1181 } else if (dir == 1) {
1182 crse(i,j,0) = Real(0.125) * (fine(ii-1,jj ,0) +
1183 fine(ii ,jj ,0) * Real(2.0) +
1184 fine(ii+1,jj ,0) +
1185 fine(ii-1,jj+1,0) +
1186 fine(ii ,jj+1,0) * Real(2.0) +
1187 fine(ii+1,jj+1,0) );
1188 } else {
1189 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)
1190 + Real(2.)*fine(ii-1,jj ,0) + Real(4.)*fine(ii ,jj ,0) + Real(2.)*fine(ii+1,jj ,0)
1191 + fine(ii-1,jj+1,0) + Real(2.)*fine(ii ,jj+1,0) + fine(ii+1,jj+1,0));
1192
1193 }
1194
1195#elif (AMREX_SPACEDIM == 1)
1196 if (dir == 0) {
1197 crse(i,0,0) = Real(0.5) * (fine(ii,0,0) + fine(ii+1,0,0));
1198 } else {
1199 crse(i,0,0) = Real(0.25) * (fine(ii-1,0,0) +
1200 fine(ii ,0,0) * Real(2.0) +
1201 fine(ii+1,0,0));
1202 }
1203#endif
1204 }
1205}
1206
1207AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1208void mlcurlcurl_bc_symmetry (int i, int j, int k, Orientation face, IndexType it,
1209 Array4<Real> const& a)
1210{
1211 int const idir = face.coordDir();
1212 int offset = face.isLow() ? 1 : -1;
1213 Real sign;
1214 if (it.cellCentered(idir)) {
1215 sign = Real(-1.0);
1216 } else {
1217 sign = Real(1.0);
1218 offset *= 2;
1219 }
1220
1221 if (idir == 0) {
1222 a(i,j,k) = sign * a(i+offset,j,k);
1223 } else if (idir == 1) {
1224 a(i,j,k) = sign * a(i,j+offset,k);
1225 } else {
1226 a(i,j,k) = sign * a(i,j,k+offset);
1227 }
1228}
1229
1230}
1231
1232#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::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, bool valid_x)
Definition AMReX_MLCurlCurl_K.H:463
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:1048
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:127
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:1208
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:1101
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