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