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 
184 struct 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 == 2)
193  amrex::ignore_unused(k);
194  return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
195  || (j == dirichlet_lo[1]) || (j == dirichlet_hi[1]);
196 #else
197  return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
198  || (j == dirichlet_lo[1]) || (j == dirichlet_hi[1])
199  || (k == dirichlet_lo[2]) || (k == dirichlet_hi[2]);
200 #endif
201  }
202 
203  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
204  bool is_dirichlet_x_edge (int, int j, int k) const
205  {
206 #if (AMREX_SPACEDIM == 2)
207  amrex::ignore_unused(k);
208  return (j == dirichlet_lo[1]) || (j == dirichlet_hi[1]);
209 #else
210  return (j == dirichlet_lo[1]) || (j == dirichlet_hi[1])
211  || (k == dirichlet_lo[2]) || (k == dirichlet_hi[2]);
212 #endif
213  }
214 
215  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
216  bool is_dirichlet_y_edge (int i, int, int k) const
217  {
218 #if (AMREX_SPACEDIM == 2)
219  amrex::ignore_unused(k);
220  return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0]);
221 #else
222  return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
223  || (k == dirichlet_lo[2]) || (k == dirichlet_hi[2]);
224 #endif
225  }
226 
227  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
228  bool is_dirichlet_z_edge (int i, int j, int) const
229  {
230  return (i == dirichlet_lo[0]) || (i == dirichlet_hi[0])
231  || (j == dirichlet_lo[1]) || (j == dirichlet_hi[1]);
232  }
233 
234  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
235  bool is_dirichlet_edge (int dim, int i, int j, int k) const
236  {
237  if (dim == 0) {
238  return is_dirichlet_x_edge(i,j,k);
239  } else if (dim == 1) {
240  return is_dirichlet_y_edge(i,j,k);
241  } else {
242  return is_dirichlet_z_edge(i,j,k);
243  }
244  }
245 };
246 
247 struct CurlCurlSymmetryInfo
248 {
249  IntVect symmetry_lo;
250  IntVect symmetry_hi;
251 
252  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
253  bool xlo_is_symmetric (int i) const
254  {
255  return i == symmetry_lo[0];
256  }
257 
258  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
259  bool xhi_is_symmetric (int i) const
260  {
261  return i == symmetry_hi[0];
262  }
263 
264  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
265  bool ylo_is_symmetric (int j) const
266  {
267  return j == symmetry_lo[1];
268  }
269 
270  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
271  bool yhi_is_symmetric (int j) const
272  {
273  return j == symmetry_hi[1];
274  }
275 
276 #if (AMREX_SPACEDIM == 3)
277  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
278  bool zlo_is_symmetric (int k) const
279  {
280  return k == symmetry_lo[2];
281  }
282 
283  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
284  bool zhi_is_symmetric (int k) const
285  {
286  return k == symmetry_hi[2];
287  }
288 #endif
289 
290  [[nodiscard]] bool is_symmetric (int dir, int side, int idx) const
291  {
292 #if (AMREX_SPACEDIM == 2)
293  if (dir == 0) {
294  return (side == 0) ? xlo_is_symmetric(idx) : xhi_is_symmetric(idx);
295  } else {
296  return (side == 0) ? ylo_is_symmetric(idx) : yhi_is_symmetric(idx);
297  }
298 #else
299  if (dir == 0) {
300  return (side == 0) ? xlo_is_symmetric(idx) : xhi_is_symmetric(idx);
301  } else if (dir == 1) {
302  return (side == 0) ? ylo_is_symmetric(idx) : yhi_is_symmetric(idx);
303  } else {
304  return (side == 0) ? zlo_is_symmetric(idx) : zhi_is_symmetric(idx);
305  }
306 #endif
307  }
308 };
309 
310 enum struct CurlCurlStateType { x, b, r }; // x, b & r=b-Ax
311 
312 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
313 void mlcurlcurl_adotx_x (int i, int j, int k, Array4<Real> const& Ax,
314  Array4<Real const> const& ex,
315  Array4<Real const> const& ey,
316  Array4<Real const> const& ez,
317  Real beta, GpuArray<Real,AMREX_SPACEDIM> const& adxinv)
318 {
319 #if (AMREX_SPACEDIM == 2)
320  amrex::ignore_unused(ez);
321  Real dyy = adxinv[1] * adxinv[1];
322  Real dxy = adxinv[0] * adxinv[1];
323  Real ccex = ex(i ,j ,k ) * dyy * Real(2.0)
324  - dyy * (ex(i ,j-1,k ) +
325  ex(i ,j+1,k ))
326  + dxy * (ey(i ,j-1,k )
327  - ey(i ,j ,k )
328  - ey(i+1,j-1,k )
329  + ey(i+1,j ,k ));
330 #else
331  Real dyy = adxinv[1] * adxinv[1];
332  Real dzz = adxinv[2] * adxinv[2];
333  Real dxy = adxinv[0] * adxinv[1];
334  Real dxz = adxinv[0] * adxinv[2];
335  Real ccex = ex(i ,j ,k ) * (dyy+dzz)*Real(2.0)
336  - dyy * (ex(i ,j-1,k ) +
337  ex(i ,j+1,k ))
338  - dzz * (ex(i ,j ,k+1) +
339  ex(i ,j ,k-1))
340  + dxy * (ey(i ,j-1,k )
341  - ey(i ,j ,k )
342  - ey(i+1,j-1,k )
343  + ey(i+1,j ,k ))
344  + dxz * (ez(i ,j ,k-1)
345  - ez(i ,j ,k )
346  - ez(i+1,j ,k-1)
347  + ez(i+1,j ,k ));
348 #endif
349  Ax(i,j,k) = ccex + beta * ex(i,j,k);
350 }
351 
352 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
353 void mlcurlcurl_adotx_y (int i, int j, int k, Array4<Real> const& Ay,
354  Array4<Real const> const& ex,
355  Array4<Real const> const& ey,
356  Array4<Real const> const& ez,
357  Real beta, GpuArray<Real,AMREX_SPACEDIM> const& adxinv)
358 {
359 #if (AMREX_SPACEDIM == 2)
360  amrex::ignore_unused(ez);
361  Real dxx = adxinv[0] * adxinv[0];
362  Real dxy = adxinv[0] * adxinv[1];
363  Real ccey = ey(i ,j ,k ) * dxx * Real(2.0)
364  - dxx * (ey(i-1,j ,k ) +
365  ey(i+1,j ,k ))
366  + dxy * (ex(i-1,j ,k )
367  - ex(i ,j ,k )
368  - ex(i-1,j+1,k )
369  + ex(i ,j+1,k ));
370 #else
371  Real dxx = adxinv[0] * adxinv[0];
372  Real dzz = adxinv[2] * adxinv[2];
373  Real dxy = adxinv[0] * adxinv[1];
374  Real dyz = adxinv[1] * adxinv[2];
375  Real ccey = ey(i ,j ,k ) * (dxx+dzz)*Real(2.0)
376  - dxx * (ey(i-1,j ,k ) +
377  ey(i+1,j ,k ))
378  - dzz * (ey(i ,j ,k-1) +
379  ey(i ,j ,k+1))
380  + dxy * (ex(i-1,j ,k )
381  - ex(i ,j ,k )
382  - ex(i-1,j+1,k )
383  + ex(i ,j+1,k ))
384  + dyz * (ez(i ,j ,k-1)
385  - ez(i ,j ,k )
386  - ez(i ,j+1,k-1)
387  + ez(i ,j+1,k ));
388 #endif
389  Ay(i,j,k) = ccey + beta * ey(i,j,k);
390 }
391 
392 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
393 void mlcurlcurl_adotx_z (int i, int j, int k, Array4<Real> const& Az,
394  Array4<Real const> const& ex,
395  Array4<Real const> const& ey,
396  Array4<Real const> const& ez,
397  Real beta, GpuArray<Real,AMREX_SPACEDIM> const& adxinv)
398 {
399 #if (AMREX_SPACEDIM == 2)
400  amrex::ignore_unused(ex,ey);
401  Real dxx = adxinv[0] * adxinv[0];
402  Real dyy = adxinv[1] * adxinv[1];
403  Real ccez = ez(i ,j ,k ) * (dxx+dyy)*Real(2.0)
404  - dxx * (ez(i-1,j ,k ) +
405  ez(i+1,j ,k ))
406  - dyy * (ez(i ,j-1,k ) +
407  ez(i ,j+1,k ));
408 #else
409  Real dxx = adxinv[0] * adxinv[0];
410  Real dyy = adxinv[1] * adxinv[1];
411  Real dxz = adxinv[0] * adxinv[2];
412  Real dyz = adxinv[1] * adxinv[2];
413  Real ccez = ez(i ,j ,k ) * (dxx+dyy)*Real(2.0)
414  - dxx * (ez(i-1,j ,k ) +
415  ez(i+1,j ,k ))
416  - dyy * (ez(i ,j-1,k ) +
417  ez(i ,j+1,k ))
418  + dxz * (ex(i-1,j ,k )
419  - ex(i ,j ,k )
420  - ex(i-1,j ,k+1)
421  + ex(i ,j ,k+1))
422  + dyz * (ey(i ,j-1,k )
423  - ey(i ,j ,k )
424  - ey(i ,j-1,k+1)
425  + ey(i ,j ,k+1));
426 #endif
427  Az(i,j,k) = ccez + beta * ez(i,j,k);
428 }
429 
430 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
431 void mlcurlcurl_gs4_lu (int i, int j, int k,
432  Array4<Real> const& ex,
433  Array4<Real> const& ey,
434  Array4<Real> const& ez,
435  Array4<Real const> const& rhsx,
436  Array4<Real const> const& rhsy,
437  Array4<Real const> const& rhsz,
438 #if (AMREX_SPACEDIM == 2)
439  Real beta,
440 #endif
441  GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
442  int color, LUSolver<AMREX_SPACEDIM*2,Real> const& lusolver,
443  CurlCurlDirichletInfo const& dinfo,
444  CurlCurlSymmetryInfo const& sinfo)
445 {
446  if (dinfo.is_dirichlet_node(i,j,k)) { return; }
447 
448  int my_color = i%2 + 2*(j%2);
449  if (k%2 != 0) {
450  my_color = 3 - my_color;
451  }
452 
453 #if (AMREX_SPACEDIM == 2)
454 
455  Real dxx = adxinv[0] * adxinv[0];
456  Real dyy = adxinv[1] * adxinv[1];
457 
458  if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) ||
459  ((my_color == 1 || my_color == 2) && (color == 1 || color == 2)))
460  {
461  Real gamma = (dxx+dyy)*Real(2.0) + beta;
462  Real ccez = - dxx * (ez(i-1,j ,k ) +
463  ez(i+1,j ,k ))
464  - dyy * (ez(i ,j-1,k ) +
465  ez(i ,j+1,k ));
466  Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez);
467  constexpr Real omega = Real(1.15);
468  ez(i,j,k) += omega/gamma * res;
469  }
470 
471  if (my_color != color) { return; }
472 
473  Real dxy = adxinv[0]*adxinv[1];
474 
475  GpuArray<Real,6> b
476  {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
477  ex(i-1,j+1,k ))
478  + dxy * ( ey(i-1,j-1,k )
479  -ey(i-1,j ,k ))),
480  rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
481  ex(i ,j+1,k ))
482  + dxy * (-ey(i+1,j-1,k )
483  +ey(i+1,j ,k ))),
484  rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
485  ey(i+1,j-1,k ))
486  + dxy * ( ex(i-1,j-1,k )
487  -ex(i ,j-1,k ))),
488  rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
489  ey(i+1,j ,k ))
490  + dxy * (-ex(i-1,j+1,k )
491  +ex(i ,j+1,k )))};
492 
493  if (sinfo.xlo_is_symmetric(i)) {
494  b[0] = -b[1];
495  } else if (sinfo.xhi_is_symmetric(i)) {
496  b[1] = -b[0];
497  }
498 
499  if (sinfo.ylo_is_symmetric(j)) {
500  b[2] = -b[3];
501  } else if (sinfo.yhi_is_symmetric(j)) {
502  b[3] = -b[2];
503  }
504 
505  GpuArray<Real,4> x;
506  lusolver(x.data(), b.data());
507  ex(i-1,j ,k ) = x[0];
508  ex(i ,j ,k ) = x[1];
509  ey(i ,j-1,k ) = x[2];
510  ey(i ,j ,k ) = x[3];
511 
512 #else
513 
514  if (my_color != color) { return; }
515 
516  Real dxx = adxinv[0]*adxinv[0];
517  Real dyy = adxinv[1]*adxinv[1];
518  Real dzz = adxinv[2]*adxinv[2];
519  Real dxy = adxinv[0]*adxinv[1];
520  Real dxz = adxinv[0]*adxinv[2];
521  Real dyz = adxinv[1]*adxinv[2];
522 
523  GpuArray<Real,6> b
524  {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
525  ex(i-1,j+1,k ))
526  - dzz * ( ex(i-1,j ,k+1) +
527  ex(i-1,j ,k-1))
528  + dxy * ( ey(i-1,j-1,k )
529  -ey(i-1,j ,k ))
530  + dxz * ( ez(i-1,j ,k-1)
531  -ez(i-1,j ,k ))),
532  rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
533  ex(i ,j+1,k ))
534  - dzz * ( ex(i ,j ,k+1) +
535  ex(i ,j ,k-1))
536  + dxy * (-ey(i+1,j-1,k )
537  +ey(i+1,j ,k ))
538  + dxz * (-ez(i+1,j ,k-1)
539  +ez(i+1,j ,k ))),
540  rhsy(i,j-1,k) - (-dxx * ( ey(i-1,j-1,k ) +
541  ey(i+1,j-1,k ))
542  - dzz * ( ey(i ,j-1,k-1) +
543  ey(i ,j-1,k+1))
544  + dxy * ( ex(i-1,j-1,k )
545  -ex(i ,j-1,k ))
546  + dyz * ( ez(i ,j-1,k-1)
547  -ez(i ,j-1,k ))),
548  rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
549  ey(i+1,j ,k ))
550  - dzz * ( ey(i ,j ,k-1) +
551  ey(i ,j ,k+1))
552  + dxy * (-ex(i-1,j+1,k )
553  +ex(i ,j+1,k ))
554  + dyz * (-ez(i ,j+1,k-1)
555  +ez(i ,j+1,k ))),
556  rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) +
557  ez(i+1,j ,k-1))
558  - dyy * ( ez(i ,j-1,k-1) +
559  ez(i ,j+1,k-1))
560  + dxz * ( ex(i-1,j ,k-1)
561  -ex(i ,j ,k-1))
562  + dyz * ( ey(i ,j-1,k-1)
563  -ey(i ,j ,k-1))),
564  rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) +
565  ez(i+1,j ,k ))
566  - dyy * ( ez(i ,j-1,k ) +
567  ez(i ,j+1,k ))
568  + dxz * (-ex(i-1,j ,k+1)
569  +ex(i ,j ,k+1))
570  + dyz * (-ey(i ,j-1,k+1)
571  +ey(i ,j ,k+1)))};
572 
573  if (sinfo.xlo_is_symmetric(i)) {
574  b[0] = -b[1];
575  } else if (sinfo.xhi_is_symmetric(i)) {
576  b[1] = -b[0];
577  }
578 
579  if (sinfo.ylo_is_symmetric(j)) {
580  b[2] = -b[3];
581  } else if (sinfo.yhi_is_symmetric(j)) {
582  b[3] = -b[2];
583  }
584 
585  if (sinfo.zlo_is_symmetric(k)) {
586  b[4] = -b[5];
587  } else if (sinfo.zhi_is_symmetric(k)) {
588  b[5] = -b[4];
589  }
590 
591  GpuArray<Real,6> x;
592  lusolver(x.data(), b.data());
593  ex(i-1,j ,k ) = x[0];
594  ex(i ,j ,k ) = x[1];
595  ey(i ,j-1,k ) = x[2];
596  ey(i ,j ,k ) = x[3];
597  ez(i ,j ,k-1) = x[4];
598  ez(i ,j ,k ) = x[5];
599 #endif
600 }
601 
602 template <bool PCG>
603 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
604 void mlcurlcurl_gs4 (int i, int j, int k,
605  Array4<Real> const& ex,
606  Array4<Real> const& ey,
607  Array4<Real> const& ez,
608  Array4<Real const> const& rhsx,
609  Array4<Real const> const& rhsy,
610  Array4<Real const> const& rhsz,
611  GpuArray<Real,AMREX_SPACEDIM> const& adxinv,
612  int color,
613  Array4<Real const> const& betax,
614  Array4<Real const> const& betay,
615  Array4<Real const> const& betaz,
616  CurlCurlDirichletInfo const& dinfo,
617  CurlCurlSymmetryInfo const& sinfo)
618 {
619  if (dinfo.is_dirichlet_node(i,j,k)) { return; }
620 
621  int my_color = i%2 + 2*(j%2);
622  if (k%2 != 0) {
623  my_color = 3 - my_color;
624  }
625 
626 #if (AMREX_SPACEDIM == 2)
627 
628  Real dxx = adxinv[0] * adxinv[0];
629  Real dyy = adxinv[1] * adxinv[1];
630 
631  if (((my_color == 0 || my_color == 3) && (color == 0 || color == 3)) ||
632  ((my_color == 1 || my_color == 2) && (color == 1 || color == 2)))
633  {
634  Real gamma = (dxx+dyy)*Real(2.0) + betaz(i,j,k);
635  Real ccez = - dxx * (ez(i-1,j ,k ) +
636  ez(i+1,j ,k ))
637  - dyy * (ez(i ,j-1,k ) +
638  ez(i ,j+1,k ));
639  Real res = rhsz(i,j,k) - (gamma*ez(i,j,k) + ccez);
640  constexpr Real omega = Real(1.15);
641  ez(i,j,k) += omega/gamma * res;
642  }
643 
644  if (my_color != color) { return; }
645 
646  Real dxy = adxinv[0]*adxinv[1];
647 
648  GpuArray<Real,6> b
649  {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
650  ex(i-1,j+1,k ))
651  + dxy * ( ey(i-1,j-1,k )
652  -ey(i-1,j ,k ))),
653  rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
654  ex(i ,j+1,k ))
655  + dxy * (-ey(i+1,j-1,k )
656  +ey(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  + dxy * ( ex(i-1,j-1,k )
660  -ex(i ,j-1,k ))),
661  rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
662  ey(i+1,j ,k ))
663  + dxy * (-ex(i-1,j+1,k )
664  +ex(i ,j+1,k )))};
665 
666  GpuArray<Real,4> beta;
667 
668  if (sinfo.xlo_is_symmetric(i)) {
669  b[0] = -b[1];
670  beta[0] = beta[1] = betax(i,j,k);
671  } else if (sinfo.xhi_is_symmetric(i)) {
672  b[1] = -b[0];
673  beta[0] = beta[1] = betax(i-1,j,k);
674  } else {
675  beta[0] = betax(i-1,j,k);
676  beta[1] = betax(i ,j,k);
677  }
678 
679  if (sinfo.ylo_is_symmetric(j)) {
680  b[2] = -b[3];
681  beta[2] = beta[3] = betay(i,j,k);
682  } else if (sinfo.yhi_is_symmetric(j)) {
683  b[3] = -b[2];
684  beta[2] = beta[3] = betay(i,j-1,k);
685  } else {
686  beta[2] = betay(i,j-1,k);
687  beta[3] = betay(i,j ,k);
688  }
689 
690  if constexpr (PCG) {
691  Real diagInv[4] = {Real(1.0) / (dyy*Real(2.0) + beta[0]),
692  Real(1.0) / (dyy*Real(2.0) + beta[1]),
693  Real(1.0) / (dxx*Real(2.0) + beta[2]),
694  Real(1.0) / (dxx*Real(2.0) + beta[3])};
695  auto precond = [&] (Real * AMREX_RESTRICT z,
696  Real const* AMREX_RESTRICT r)
697  {
698  for (int m = 0; m < 4; ++m) { z[m] = r[m] * diagInv[m]; }
699  };
700  auto mat = [&] (Real * AMREX_RESTRICT Av,
701  Real const* AMREX_RESTRICT v)
702  {
703  Av[0] = (dyy*Real(2.0) + beta[0]) * v[0] - dxy * v[2] + dxy * v[3];
704  Av[1] = (dyy*Real(2.0) + beta[1]) * v[1] + dxy * v[2] - dxy * v[3];
705  Av[2] = -dxy * v[0] + dxy * v[1] + (dxx*Real(2.0) + beta[2]) * v[2];
706  Av[3] = dxy * v[0] - dxy * v[1] + (dxx*Real(2.0) + beta[3]) * v[3];
707  };
708  Real sol[4] = {0, 0, 0, 0};
709  pcg_solve<4>(sol, b.data(), mat, precond, 8, Real(1.e-8));
710  ex(i-1,j ,k ) = sol[0];
711  ex(i ,j ,k ) = sol[1];
712  ey(i ,j-1,k ) = sol[2];
713  ey(i ,j ,k ) = sol[3];
714  } else {
715  LUSolver<4,Real> lusolver
716  ({dyy*Real(2.0) + beta[0],
717  Real(0.0),
718  -dxy,
719  dxy,
720  //
721  Real(0.0),
722  dyy*Real(2.0) + beta[1],
723  dxy,
724  -dxy,
725  //
726  -dxy,
727  dxy,
728  dxx*Real(2.0) + beta[2],
729  Real(0.0),
730  //
731  dxy,
732  -dxy,
733  Real(0.0),
734  dxx*Real(2.0) + beta[3]});
735  lusolver(beta.data(), b.data());
736  ex(i-1,j ,k ) = beta[0];
737  ex(i ,j ,k ) = beta[1];
738  ey(i ,j-1,k ) = beta[2];
739  ey(i ,j ,k ) = beta[3];
740  }
741 
742 #else
743 
744  if (my_color != color) { return; }
745 
746  Real dxx = adxinv[0]*adxinv[0];
747  Real dyy = adxinv[1]*adxinv[1];
748  Real dzz = adxinv[2]*adxinv[2];
749  Real dxy = adxinv[0]*adxinv[1];
750  Real dxz = adxinv[0]*adxinv[2];
751  Real dyz = adxinv[1]*adxinv[2];
752 
753  GpuArray<Real,6> b
754  {rhsx(i-1,j,k) - (-dyy * ( ex(i-1,j-1,k ) +
755  ex(i-1,j+1,k ))
756  - dzz * ( ex(i-1,j ,k+1) +
757  ex(i-1,j ,k-1))
758  + dxy * ( ey(i-1,j-1,k )
759  -ey(i-1,j ,k ))
760  + dxz * ( ez(i-1,j ,k-1)
761  -ez(i-1,j ,k ))),
762  rhsx(i ,j,k) - (-dyy * ( ex(i ,j-1,k ) +
763  ex(i ,j+1,k ))
764  - dzz * ( ex(i ,j ,k+1) +
765  ex(i ,j ,k-1))
766  + dxy * (-ey(i+1,j-1,k )
767  +ey(i+1,j ,k ))
768  + dxz * (-ez(i+1,j ,k-1)
769  +ez(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  - dzz * ( ey(i ,j-1,k-1) +
773  ey(i ,j-1,k+1))
774  + dxy * ( ex(i-1,j-1,k )
775  -ex(i ,j-1,k ))
776  + dyz * ( ez(i ,j-1,k-1)
777  -ez(i ,j-1,k ))),
778  rhsy(i,j ,k) - (-dxx * ( ey(i-1,j ,k ) +
779  ey(i+1,j ,k ))
780  - dzz * ( ey(i ,j ,k-1) +
781  ey(i ,j ,k+1))
782  + dxy * (-ex(i-1,j+1,k )
783  +ex(i ,j+1,k ))
784  + dyz * (-ez(i ,j+1,k-1)
785  +ez(i ,j+1,k ))),
786  rhsz(i,j,k-1) - (-dxx * ( ez(i-1,j ,k-1) +
787  ez(i+1,j ,k-1))
788  - dyy * ( ez(i ,j-1,k-1) +
789  ez(i ,j+1,k-1))
790  + dxz * ( ex(i-1,j ,k-1)
791  -ex(i ,j ,k-1))
792  + dyz * ( ey(i ,j-1,k-1)
793  -ey(i ,j ,k-1))),
794  rhsz(i,j,k ) - (-dxx * ( ez(i-1,j ,k ) +
795  ez(i+1,j ,k ))
796  - dyy * ( ez(i ,j-1,k ) +
797  ez(i ,j+1,k ))
798  + dxz * (-ex(i-1,j ,k+1)
799  +ex(i ,j ,k+1))
800  + dyz * (-ey(i ,j-1,k+1)
801  +ey(i ,j ,k+1)))};
802 
803  GpuArray<Real,6> beta;
804 
805  if (sinfo.xlo_is_symmetric(i)) {
806  b[0] = -b[1];
807  beta[0] = beta[1] = betax(i,j,k);
808  } else if (sinfo.xhi_is_symmetric(i)) {
809  b[1] = -b[0];
810  beta[0] = beta[1] = betax(i-1,j,k);
811  } else {
812  beta[0] = betax(i-1,j,k);
813  beta[1] = betax(i ,j,k);
814  }
815 
816  if (sinfo.ylo_is_symmetric(j)) {
817  b[2] = -b[3];
818  beta[2] = beta[3] = betay(i,j,k);
819  } else if (sinfo.yhi_is_symmetric(j)) {
820  b[3] = -b[2];
821  beta[2] = beta[3] = betay(i,j-1,k);
822  } else {
823  beta[2] = betay(i,j-1,k);
824  beta[3] = betay(i,j ,k);
825  }
826 
827  if (sinfo.zlo_is_symmetric(k)) {
828  b[4] = -b[5];
829  beta[4] = beta[5] = betaz(i,j,k);
830  } else if (sinfo.zhi_is_symmetric(k)) {
831  b[5] = -b[4];
832  beta[4] = beta[5] = betaz(i,j,k-1);
833  } else {
834  beta[4] = betaz(i,j,k-1);
835  beta[5] = betaz(i,j,k );
836  }
837 
838  if constexpr (PCG) {
839  Real diagInv[6] = {Real(1.0) / ((dyy+dzz)*Real(2.0) + beta[0]),
840  Real(1.0) / ((dyy+dzz)*Real(2.0) + beta[1]),
841  Real(1.0) / ((dxx+dzz)*Real(2.0) + beta[2]),
842  Real(1.0) / ((dxx+dzz)*Real(2.0) + beta[3]),
843  Real(1.0) / ((dxx+dyy)*Real(2.0) + beta[4]),
844  Real(1.0) / ((dxx+dyy)*Real(2.0) + beta[5])};
845  auto precond = [&] (Real * AMREX_RESTRICT z,
846  Real const* AMREX_RESTRICT r)
847  {
848  for (int m = 0; m < 6; ++m) { z[m] = r[m] * diagInv[m]; }
849  };
850  auto mat = [&] (Real * AMREX_RESTRICT Av,
851  Real const* AMREX_RESTRICT v)
852  {
853  Av[0] = ((dyy+dzz)*Real(2.0) + beta[0]) * v[0] - dxy * v[2]
854  + dxy * v[3] - dxz * v[4] + dxz * v[5];
855  Av[1] = ((dyy+dzz)*Real(2.0) + beta[1]) * v[1] + dxy * v[2]
856  - dxy * v[3] + dxz * v[4] - dxz * v[5];
857  Av[2] = -dxy * v[0] + dxy * v[1] + ((dxx+dzz)*Real(2.0) + beta[2]) * v[2]
858  - dyz * v[4] + dyz * v[5];
859  Av[3] = dxy * v[0] - dxy * v[1] + ((dxx+dzz)*Real(2.0) + beta[3]) * v[3]
860  + dyz * v[4] - dyz * v[5];
861  Av[4] = -dxz * v[0] + dxz * v[1] - dyz * v[2] + dyz * v[3]
862  + ((dxx+dyy)*Real(2.0) + beta[4]) * v[4];
863  Av[5] = dxz * v[0] - dxz * v[1] + dyz * v[2] - dyz * v[3]
864  + ((dxx+dyy)*Real(2.0) + beta[5]) * v[5];
865  };
866  Real sol[6] = {0, 0, 0, 0, 0, 0};
867  pcg_solve<6>(sol, b.data(), mat, precond, 8, Real(1.e-8));
868  ex(i-1,j ,k ) = sol[0];
869  ex(i ,j ,k ) = sol[1];
870  ey(i ,j-1,k ) = sol[2];
871  ey(i ,j ,k ) = sol[3];
872  ez(i ,j ,k-1) = sol[4];
873  ez(i ,j ,k ) = sol[5];
874  } else {
875  LUSolver<6,Real> lusolver
876  ({(dyy+dzz)*Real(2.0) + beta[0],
877  Real(0.0),
878  -dxy,
879  dxy,
880  -dxz,
881  dxz,
882  //
883  Real(0.0),
884  (dyy+dzz)*Real(2.0) + beta[1],
885  dxy,
886  -dxy,
887  dxz,
888  -dxz,
889  //
890  -dxy,
891  dxy,
892  (dxx+dzz)*Real(2.0) + beta[2],
893  Real(0.0),
894  -dyz,
895  dyz,
896  //
897  dxy,
898  -dxy,
899  Real(0.0),
900  (dxx+dzz)*Real(2.0) + beta[3],
901  dyz,
902  -dyz,
903  //
904  -dxz,
905  dxz,
906  -dyz,
907  dyz,
908  (dxx+dyy)*Real(2.0) + beta[4],
909  Real(0.0),
910  //
911  dxz,
912  -dxz,
913  dyz,
914  -dyz,
915  Real(0.0),
916  (dxx+dyy)*Real(2.0) + beta[5]});
917  lusolver(beta.data(), b.data());
918  ex(i-1,j ,k ) = beta[0];
919  ex(i ,j ,k ) = beta[1];
920  ey(i ,j-1,k ) = beta[2];
921  ey(i ,j ,k ) = beta[3];
922  ez(i ,j ,k-1) = beta[4];
923  ez(i ,j ,k ) = beta[5];
924  }
925 #endif
926 }
927 
928 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
929 void mlcurlcurl_interpadd (int dir, int i, int j, int k,
930  Array4<Real> const& fine,
931  Array4<Real const> const& crse)
932 {
933  int ic = amrex::coarsen(i,2);
934  int jc = amrex::coarsen(j,2);
935  int kc = amrex::coarsen(k,2);
936  if (dir == 0) {
937  bool j_is_odd = (jc*2 != j);
938  bool k_is_odd = (kc*2 != k);
939  if (j_is_odd && k_is_odd) {
940  fine(i,j,k) += Real(0.25) *
941  (crse(ic,jc,kc ) + crse(ic,jc+1,kc ) +
942  crse(ic,jc,kc+1) + crse(ic,jc+1,kc+1));
943  } else if (j_is_odd) {
944  fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc));
945  } else if (k_is_odd) {
946  fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1));
947  } else {
948  fine(i,j,k) += crse(ic,jc,kc);
949  }
950  } else if (dir == 1) {
951  bool i_is_odd = (ic*2 != i);
952  bool k_is_odd = (kc*2 != k);
953  if (i_is_odd && k_is_odd) {
954  fine(i,j,k) += Real(0.25) *
955  (crse(ic ,jc,kc ) + crse(ic+1,jc,kc ) +
956  crse(ic ,jc,kc+1) + crse(ic+1,jc,kc+1));
957  } else if (i_is_odd) {
958  fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc));
959  } else if (k_is_odd) {
960  fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc,kc+1));
961  } else {
962  fine(i,j,k) += crse(ic,jc,kc);
963  }
964  } else {
965  bool i_is_odd = (ic*2 != i);
966  bool j_is_odd = (jc*2 != j);
967  if (i_is_odd && j_is_odd) {
968  fine(i,j,k) += Real(0.25) *
969  (crse(ic ,jc ,kc) + crse(ic+1,jc ,kc) +
970  crse(ic ,jc+1,kc) + crse(ic+1,jc+1,kc));
971  } else if (i_is_odd) {
972  fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic+1,jc,kc));
973  } else if (j_is_odd) {
974  fine(i,j,k) += Real(0.5) * (crse(ic,jc,kc) + crse(ic,jc+1,kc));
975  } else {
976  fine(i,j,k) += crse(ic,jc,kc);
977  }
978  }
979 }
980 
981 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
982 void mlcurlcurl_restriction (int dir, int i, int j, int k,
983  Array4<Real> const& crse,
984  Array4<Real const> const& fine,
985  CurlCurlDirichletInfo const& dinfo)
986 {
987  int ii = i*2;
988  int jj = j*2;
989  int kk = k*2;
990  if (dinfo.is_dirichlet_edge(dir,ii,jj,kk)) {
991  crse(i,j,k) = Real(0.0);
992  }
993  else
994  {
995 #if (AMREX_SPACEDIM == 3)
996  if (dir == 0) {
997  crse(i,j,k) = Real(1./32.) * (fine(ii ,jj-1,kk-1) +
998  fine(ii ,jj ,kk-1) * Real(2.0) +
999  fine(ii ,jj+1,kk-1) +
1000  fine(ii ,jj-1,kk ) * Real(2.0) +
1001  fine(ii ,jj ,kk ) * Real(4.0) +
1002  fine(ii ,jj+1,kk ) * Real(2.0) +
1003  fine(ii ,jj-1,kk+1) +
1004  fine(ii ,jj ,kk+1) * Real(2.0) +
1005  fine(ii ,jj+1,kk+1) +
1006  fine(ii+1,jj-1,kk-1) +
1007  fine(ii+1,jj ,kk-1) * Real(2.0) +
1008  fine(ii+1,jj+1,kk-1) +
1009  fine(ii+1,jj-1,kk ) * Real(2.0) +
1010  fine(ii+1,jj ,kk ) * Real(4.0) +
1011  fine(ii+1,jj+1,kk ) * Real(2.0) +
1012  fine(ii+1,jj-1,kk+1) +
1013  fine(ii+1,jj ,kk+1) * Real(2.0) +
1014  fine(ii+1,jj+1,kk+1) );
1015  } else if (dir == 1) {
1016  crse(i,j,k) = Real(1./32.) * (fine(ii-1,jj ,kk-1) +
1017  fine(ii ,jj ,kk-1) * Real(2.0) +
1018  fine(ii+1,jj ,kk-1) +
1019  fine(ii-1,jj ,kk ) * Real(2.0) +
1020  fine(ii ,jj ,kk ) * Real(4.0) +
1021  fine(ii+1,jj ,kk ) * Real(2.0) +
1022  fine(ii-1,jj ,kk+1) +
1023  fine(ii ,jj ,kk+1) * Real(2.0) +
1024  fine(ii+1,jj ,kk+1) +
1025  fine(ii-1,jj+1,kk-1) +
1026  fine(ii ,jj+1,kk-1) * Real(2.0) +
1027  fine(ii+1,jj+1,kk-1) +
1028  fine(ii-1,jj+1,kk ) * Real(2.0) +
1029  fine(ii ,jj+1,kk ) * Real(4.0) +
1030  fine(ii+1,jj+1,kk ) * Real(2.0) +
1031  fine(ii-1,jj+1,kk+1) +
1032  fine(ii ,jj+1,kk+1) * Real(2.0) +
1033  fine(ii+1,jj+1,kk+1) );
1034  } else {
1035  crse(i,j,k) = Real(1./32.) * (fine(ii-1,jj-1,kk ) +
1036  fine(ii ,jj-1,kk ) * Real(2.0) +
1037  fine(ii+1,jj-1,kk ) +
1038  fine(ii-1,jj ,kk ) * Real(2.0) +
1039  fine(ii ,jj ,kk ) * Real(4.0) +
1040  fine(ii+1,jj ,kk ) * Real(2.0) +
1041  fine(ii-1,jj+1,kk ) +
1042  fine(ii ,jj+1,kk ) * Real(2.0) +
1043  fine(ii+1,jj+1,kk ) +
1044  fine(ii-1,jj-1,kk+1) +
1045  fine(ii ,jj-1,kk+1) * Real(2.0) +
1046  fine(ii+1,jj-1,kk+1) +
1047  fine(ii-1,jj ,kk+1) * Real(2.0) +
1048  fine(ii ,jj ,kk+1) * Real(4.0) +
1049  fine(ii+1,jj ,kk+1) * Real(2.0) +
1050  fine(ii-1,jj+1,kk+1) +
1051  fine(ii ,jj+1,kk+1) * Real(2.0) +
1052  fine(ii+1,jj+1,kk+1) );
1053  }
1054 #else
1055  if (dir == 0) {
1056  crse(i,j,0) = Real(0.125) * (fine(ii ,jj-1,0) +
1057  fine(ii ,jj ,0) * Real(2.0) +
1058  fine(ii ,jj+1,0) +
1059  fine(ii+1,jj-1,0) +
1060  fine(ii+1,jj ,0) * Real(2.0) +
1061  fine(ii+1,jj+1,0) );
1062  } else if (dir == 1) {
1063  crse(i,j,0) = Real(0.125) * (fine(ii-1,jj ,0) +
1064  fine(ii ,jj ,0) * Real(2.0) +
1065  fine(ii+1,jj ,0) +
1066  fine(ii-1,jj+1,0) +
1067  fine(ii ,jj+1,0) * Real(2.0) +
1068  fine(ii+1,jj+1,0) );
1069  } else {
1070  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)
1071  + Real(2.)*fine(ii-1,jj ,0) + Real(4.)*fine(ii ,jj ,0) + Real(2.)*fine(ii+1,jj ,0)
1072  + fine(ii-1,jj+1,0) + Real(2.)*fine(ii ,jj+1,0) + fine(ii+1,jj+1,0));
1073 
1074  }
1075 #endif
1076  }
1077 }
1078 
1079 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
1080 void mlcurlcurl_bc_symmetry (int i, int j, int k, Orientation face, IndexType it,
1081  Array4<Real> const& a)
1082 {
1083  int const idir = face.coordDir();
1084  int offset = face.isLow() ? 1 : -1;
1085  Real sign;
1086  if (it.cellCentered(idir)) {
1087  sign = Real(-1.0);
1088  } else {
1089  sign = Real(1.0);
1090  offset *= 2;
1091  }
1092 
1093  if (idir == 0) {
1094  a(i,j,k) = sign * a(i+offset,j,k);
1095  } else if (idir == 1) {
1096  a(i,j,k) = sign * a(i,j+offset,k);
1097  } else {
1098  a(i,j,k) = sign * a(i,j,k+offset);
1099  }
1100 }
1101 
1102 }
1103 
1104 #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::LUSolver
Definition: AMReX_LUSolver.H:16
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:310
amrex::mlcurlcurl_gs4
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_gs4(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, GpuArray< Real, AMREX_SPACEDIM > const &adxinv, int color, Array4< Real const > const &betax, Array4< Real const > const &betay, Array4< Real const > const &betaz, CurlCurlDirichletInfo const &dinfo, CurlCurlSymmetryInfo const &sinfo)
Definition: AMReX_MLCurlCurl_K.H:604
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:929
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:353
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:1080
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:393
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:313
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_gs4_lu
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void mlcurlcurl_gs4_lu(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, GpuArray< Real, AMREX_SPACEDIM > const &adxinv, int color, LUSolver< AMREX_SPACEDIM *2, Real > const &lusolver, CurlCurlDirichletInfo const &dinfo, CurlCurlSymmetryInfo const &sinfo)
Definition: AMReX_MLCurlCurl_K.H:431
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:982
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:235
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:228
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:216
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:204
amrex::CurlCurlSymmetryInfo
Definition: AMReX_MLCurlCurl_K.H:248
amrex::CurlCurlSymmetryInfo::xlo_is_symmetric
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool xlo_is_symmetric(int i) const
Definition: AMReX_MLCurlCurl_K.H:253
amrex::CurlCurlSymmetryInfo::ylo_is_symmetric
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool ylo_is_symmetric(int j) const
Definition: AMReX_MLCurlCurl_K.H:265
amrex::CurlCurlSymmetryInfo::xhi_is_symmetric
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool xhi_is_symmetric(int i) const
Definition: AMReX_MLCurlCurl_K.H:259
amrex::CurlCurlSymmetryInfo::symmetry_hi
IntVect symmetry_hi
Definition: AMReX_MLCurlCurl_K.H:250
amrex::CurlCurlSymmetryInfo::symmetry_lo
IntVect symmetry_lo
Definition: AMReX_MLCurlCurl_K.H:249
amrex::CurlCurlSymmetryInfo::yhi_is_symmetric
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool yhi_is_symmetric(int j) const
Definition: AMReX_MLCurlCurl_K.H:271
amrex::CurlCurlSymmetryInfo::is_symmetric
bool is_symmetric(int dir, int side, int idx) const
Definition: AMReX_MLCurlCurl_K.H:290
amrex::GpuArray::data
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * data() const noexcept
Definition: AMReX_Array.H:56