Block-Structured AMR Software Framework
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 
9 namespace 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 
185 {
188 
190  bool is_dirichlet_node (int i, int j, int k) const
191  {
192 #if (AMREX_SPACEDIM == 1)
194  return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0]);
195 #elif (AMREX_SPACEDIM == 2)
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 
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)
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 
222  bool is_dirichlet_y_edge (int i, int, int k) const
223  {
224 #if (AMREX_SPACEDIM < 3)
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 
234  bool is_dirichlet_z_edge (int i, int j, int) const
235  {
236 #if (AMREX_SPACEDIM == 1)
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 
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 
259 {
262 
264  bool xlo_is_symmetric (int i) const
265  {
266  return i == symmetry_lo[0];
267  }
268 
270  bool xhi_is_symmetric (int i) const
271  {
272  return i == symmetry_hi[0];
273  }
274 
275 #if (AMREX_SPACEDIM > 1)
277  bool ylo_is_symmetric (int j) const
278  {
279  return j == symmetry_lo[1];
280  }
281 
283  bool yhi_is_symmetric (int j) const
284  {
285  return j == symmetry_hi[1];
286  }
287 #endif
288 
289 #if (AMREX_SPACEDIM == 3)
291  bool zlo_is_symmetric (int k) const
292  {
293  return k == symmetry_lo[2];
294  }
295 
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)
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 
326 enum struct CurlCurlStateType { x, b, r }; // x, b & r=b-Ax
327 
329 void 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)
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 
372 void 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)
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 
418 void 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)
463 void 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 
501 void 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 
544 void 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 
715 template <bool PCG>
717 void 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 
1044 void 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 
1097 void 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 
1204 void 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
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_RESTRICT
Definition: AMReX_Extension.H:37
#define AMREX_GPU_DEVICE
Definition: AMReX_GpuQualifiers.H:18
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
int idir
Definition: AMReX_HypreMLABecLap.cpp:1093
Array4< int const > offset
Definition: AMReX_HypreMLABecLap.cpp:1089
Array4< Real > fine
Definition: AMReX_InterpFaceRegister.cpp:90
Array4< Real const > crse
Definition: AMReX_InterpFaceRegister.cpp:92
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
Encapsulation of the Orientation of the Faces of a Box.
Definition: AMReX_Orientation.H:29
AMREX_GPU_HOST_DEVICE int coordDir() const noexcept
Returns the coordinate direction.
Definition: AMReX_Orientation.H:83
AMREX_GPU_HOST_DEVICE bool isLow() const noexcept
Returns true if Orientation is low.
Definition: AMReX_Orientation.H:89
Definition: AMReX_Amr.cpp:49
CurlCurlStateType
Definition: AMReX_MLCurlCurl_K.H:326
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_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_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_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition: AMReX.H:111
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_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_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_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_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
Definition: AMReX_Array4.H:61
Definition: AMReX_MLCurlCurl_K.H:185
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_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_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_z_edge(int i, int j, int) const
Definition: AMReX_MLCurlCurl_K.H:234
IntVect dirichlet_hi
Definition: AMReX_MLCurlCurl_K.H:187
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_y_edge(int i, int, int k) const
Definition: AMReX_MLCurlCurl_K.H:222
IntVect dirichlet_lo
Definition: AMReX_MLCurlCurl_K.H:186
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool is_dirichlet_x_edge(int, int j, int k) const
Definition: AMReX_MLCurlCurl_K.H:207
Definition: AMReX_MLCurlCurl_K.H:259
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool xlo_is_symmetric(int i) const
Definition: AMReX_MLCurlCurl_K.H:264
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool xhi_is_symmetric(int i) const
Definition: AMReX_MLCurlCurl_K.H:270
IntVect symmetry_hi
Definition: AMReX_MLCurlCurl_K.H:261
IntVect symmetry_lo
Definition: AMReX_MLCurlCurl_K.H:260
bool is_symmetric(int dir, int side, int idx) const
Definition: AMReX_MLCurlCurl_K.H:303