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