WarpX
CartesianCKCAlgorithm.H
Go to the documentation of this file.
1 /* Copyright 2020 Remi Lehe
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_
9 #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_
10 
11 #include "Utils/WarpXConst.H"
12 
13 #include <AMReX.H>
14 #include <AMReX_Array4.H>
15 #include <AMReX_Gpu.H>
16 #include <AMReX_REAL.H>
17 
18 #include <algorithm>
19 #include <array>
20 
21 
27 
29  std::array<amrex::Real,3>& cell_size,
30  amrex::Vector<amrex::Real>& stencil_coefs_x,
31  amrex::Vector<amrex::Real>& stencil_coefs_y,
32  amrex::Vector<amrex::Real>& stencil_coefs_z ) {
33 
34  using namespace amrex;
35 
36  // Compute Cole-Karkkainen-Cowan coefficients according
37  // to Cowan - PRST-AB 16, 041303 (2013)
38  Real const inv_dx = 1._rt/cell_size[0];
39  Real const inv_dy = 1._rt/cell_size[1];
40  Real const inv_dz = 1._rt/cell_size[2];
41 #if defined WARPX_DIM_3D
42  Real const delta = std::max( { inv_dx,inv_dy,inv_dz } );
43  Real const rx = (inv_dx/delta)*(inv_dx/delta);
44  Real const ry = (inv_dy/delta)*(inv_dy/delta);
45  Real const rz = (inv_dz/delta)*(inv_dz/delta);
46  Real const beta = 0.125_rt*(1._rt - rx*ry*rz/(ry*rz + rz*rx + rx*ry));
47  Real const betaxy = ry*beta*inv_dx;
48  Real const betaxz = rz*beta*inv_dx;
49  Real const betayx = rx*beta*inv_dy;
50  Real const betayz = rz*beta*inv_dy;
51  Real const betazx = rx*beta*inv_dz;
52  Real const betazy = ry*beta*inv_dz;
53  Real const inv_r_fac = (1._rt/(ry*rz + rz*rx + rx*ry));
54  Real const gammax = ry*rz*(0.0625_rt - 0.125_rt*ry*rz*inv_r_fac);
55  Real const gammay = rx*rz*(0.0625_rt - 0.125_rt*rx*rz*inv_r_fac);
56  Real const gammaz = rx*ry*(0.0625_rt - 0.125_rt*rx*ry*inv_r_fac);
57  Real const alphax = (1._rt - 2._rt*ry*beta - 2._rt*rz*beta - 4._rt*gammax)*inv_dx;
58  Real const alphay = (1._rt - 2._rt*rx*beta - 2._rt*rz*beta - 4._rt*gammay)*inv_dy;
59  Real const alphaz = (1._rt - 2._rt*rx*beta - 2._rt*ry*beta - 4._rt*gammaz)*inv_dz;
60 #elif defined WARPX_DIM_XZ
61  Real const delta = std::max(inv_dx,inv_dz);
62  Real const rx = (inv_dx/delta)*(inv_dx/delta);
63  Real const rz = (inv_dz/delta)*(inv_dz/delta);
64  constexpr Real beta = 0.125_rt;
65  Real const betaxz = beta*rz*inv_dx;
66  Real const betazx = beta*rx*inv_dz;
67  Real const alphax = (1._rt - 2._rt*rz*beta)*inv_dx;
68  Real const alphaz = (1._rt - 2._rt*rx*beta)*inv_dz;
69  // Other coefficients are 0 in 2D Cartesian
70  // (and will actually not be used in the stencil)
71  constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt;
72  constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt;
73  constexpr Real alphay=0._rt;
74 #elif defined WARPX_DIM_1D_Z
75  Real const alphaz = inv_dz;
76  // Other coefficients are 0 in 1D Cartesian
77  // (and will actually not be used in the stencil)
78  constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt;
79  constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt, betaxz=0._rt, betazx=0._rt;
80  constexpr Real alphay=0._rt, alphax=0._rt;
81 #endif
82 
83  // Store the coefficients in array `stencil_coefs`, in prescribed order
84  stencil_coefs_x.resize(6);
85  stencil_coefs_x[0] = inv_dx;
86  stencil_coefs_x[1] = alphax;
87  stencil_coefs_x[2] = betaxy;
88  stencil_coefs_x[3] = betaxz;
89  stencil_coefs_x[4] = gammax*inv_dx;
90  stencil_coefs_y.resize(6);
91  stencil_coefs_y[0] = inv_dy;
92  stencil_coefs_y[1] = alphay;
93  stencil_coefs_y[2] = betayz;
94  stencil_coefs_y[3] = betayx;
95  stencil_coefs_y[4] = gammay*inv_dy;
96  stencil_coefs_z.resize(6);
97  stencil_coefs_z[0] = inv_dz;
98  stencil_coefs_z[1] = alphaz;
99  stencil_coefs_z[2] = betazx;
100  stencil_coefs_z[3] = betazy;
101  stencil_coefs_z[4] = gammaz*inv_dz;
102  }
103 
107  static amrex::Real ComputeMaxDt ( amrex::Real const * const dx ) {
108 #if (defined WARPX_DIM_1D_Z)
109  amrex::Real const delta_t = dx[0]/PhysConst::c;
110 #elif (defined WARPX_DIM_XZ)
111  // - In Cartesian 2D geometry: determined by the minimum cell size in all direction
112  amrex::Real const delta_t = std::min( dx[0], dx[1] )/PhysConst::c;
113 #else
114  // - In Cartesian 3D geometry: determined by the minimum cell size in all direction
115  amrex::Real const delta_t = std::min( dx[0], std::min( dx[1], dx[2] ) ) / PhysConst::c;
116 #endif
117  return delta_t;
118  }
119 
124  // The ckc solver requires one guard cell in each dimension
125  return amrex::IntVect{AMREX_D_DECL(1,1,1)};
126  }
127 
131  static amrex::Real UpwardDx (
133  amrex::Real const * const coefs_x, int const /*n_coefs_x*/,
134  int const i, int const j, int const k, int const ncomp=0 ) {
135 
136  using namespace amrex;
137 #if (defined WARPX_DIM_3D || WARPX_DIM_XZ)
138  amrex::Real const alphax = coefs_x[1];
139  amrex::Real const betaxz = coefs_x[3];
140 #endif
141 #if defined WARPX_DIM_3D
142  amrex::Real const betaxy = coefs_x[2];
143  amrex::Real const gammax = coefs_x[4];
144 #endif
145 
146 #if defined WARPX_DIM_3D
147  return alphax * (F(i+1,j ,k ,ncomp) - F(i, j, k ,ncomp))
148  + betaxy * (F(i+1,j+1,k ,ncomp) - F(i ,j+1,k ,ncomp)
149  + F(i+1,j-1,k ,ncomp) - F(i ,j-1,k ,ncomp))
150  + betaxz * (F(i+1,j ,k+1,ncomp) - F(i ,j ,k+1,ncomp)
151  + F(i+1,j ,k-1,ncomp) - F(i ,j ,k-1,ncomp))
152  + gammax * (F(i+1,j+1,k+1,ncomp) - F(i ,j+1,k+1,ncomp)
153  + F(i+1,j-1,k+1,ncomp) - F(i ,j-1,k+1,ncomp)
154  + F(i+1,j+1,k-1,ncomp) - F(i ,j+1,k-1,ncomp)
155  + F(i+1,j-1,k-1,ncomp) - F(i ,j-1,k-1,ncomp));
156 #elif (defined WARPX_DIM_XZ)
157  return alphax * (F(i+1,j ,k ,ncomp) - F(i, j, k ,ncomp))
158  + betaxz * (F(i+1,j+1,k ,ncomp) - F(i ,j+1,k ,ncomp)
159  + F(i+1,j-1,k ,ncomp) - F(i ,j-1,k ,ncomp));
160 #elif (defined WARPX_DIM_1D_Z)
161  amrex::ignore_unused(F, coefs_x, i, j, k, ncomp);
162  return 0._rt; // 1D Cartesian: derivative along x is 0
163 #endif
164  }
165 
168  template< typename T_Field>
170  static amrex::Real DownwardDx (
171  T_Field const& F,
172  amrex::Real const * const coefs_x, int const /*n_coefs_x*/,
173  int const i, int const j, int const k, int const ncomp=0 ) {
174 
175  using namespace amrex;
176 #if (defined WARPX_DIM_1D_Z)
177  amrex::ignore_unused(F, coefs_x, i, j, k, ncomp);
178  return 0._rt; // 1D Cartesian: derivative along x is 0
179 #else
180  amrex::Real const inv_dx = coefs_x[0];
181  return inv_dx*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) );
182 #endif
183  }
184 
188  static amrex::Real UpwardDy (
190  amrex::Real const * const coefs_y, int const n_coefs_y,
191  int const i, int const j, int const k, int const ncomp=0 ) {
192 
193  using namespace amrex;
194 #if defined WARPX_DIM_3D
195  Real const alphay = coefs_y[1];
196  Real const betayz = coefs_y[2];
197  Real const betayx = coefs_y[3];
198  Real const gammay = coefs_y[4];
199  amrex::ignore_unused(n_coefs_y);
200  return alphay * (F(i ,j+1,k ,ncomp) - F(i ,j ,k ,ncomp))
201  + betayx * (F(i+1,j+1,k ,ncomp) - F(i+1,j ,k ,ncomp)
202  + F(i-1,j+1,k ,ncomp) - F(i-1,j ,k ,ncomp))
203  + betayz * (F(i ,j+1,k+1,ncomp) - F(i ,j ,k+1,ncomp)
204  + F(i ,j+1,k-1,ncomp) - F(i ,j ,k-1,ncomp))
205  + gammay * (F(i+1,j+1,k+1,ncomp) - F(i+1,j ,k+1,ncomp)
206  + F(i-1,j+1,k+1,ncomp) - F(i-1,j ,k+1,ncomp)
207  + F(i+1,j+1,k-1,ncomp) - F(i+1,j ,k-1,ncomp)
208  + F(i-1,j+1,k-1,ncomp) - F(i-1,j ,k-1,ncomp));
209 #elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z)
210  amrex::ignore_unused(F, coefs_y, n_coefs_y,
211  i, j, k, ncomp);
212  return 0._rt; // 1D and 2D Cartesian: derivative along y is 0
213 #else
214  amrex::ignore_unused(F, coefs_y, n_coefs_y,
215  i, j, k, ncomp);
216 #endif
217  }
218 
221  template< typename T_Field>
223  static amrex::Real DownwardDy (
224  T_Field const& F,
225  amrex::Real const * const coefs_y, int const n_coefs_y,
226  int const i, int const j, int const k, int const ncomp=0 ) {
227 
228  using namespace amrex;
229 #if defined WARPX_DIM_3D
230  amrex::Real const inv_dy = coefs_y[0];
231  amrex::ignore_unused(n_coefs_y);
232  return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) );
233 #elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z)
234  amrex::ignore_unused(F, coefs_y, n_coefs_y,
235  i, j, k, ncomp);
236  return 0._rt; // 2D Cartesian: derivative along y is 0
237 #else
238  amrex::ignore_unused(F, coefs_y, n_coefs_y,
239  i, j, k, ncomp);
240 #endif
241  }
242 
246  static amrex::Real UpwardDz (
248  amrex::Real const * const coefs_z, int const n_coefs_z,
249  int const i, int const j, int const k, int const ncomp=0 ) {
250 
251  using namespace amrex;
252 
253  amrex::ignore_unused(n_coefs_z);
254 
255  Real const alphaz = coefs_z[1];
256 #if (defined WARPX_DIM_3D || WARPX_DIM_XZ)
257  Real const betazx = coefs_z[2];
258 #endif
259 #if defined WARPX_DIM_3D
260  Real const betazy = coefs_z[3];
261  Real const gammaz = coefs_z[4];
262 #endif
263 #if defined WARPX_DIM_3D
264  return alphaz * (F(i ,j ,k+1,ncomp) - F(i ,j ,k ,ncomp))
265  + betazx * (F(i+1,j ,k+1,ncomp) - F(i+1,j ,k ,ncomp)
266  + F(i-1,j ,k+1,ncomp) - F(i-1,j ,k ,ncomp))
267  + betazy * (F(i ,j+1,k+1,ncomp) - F(i ,j+1,k ,ncomp)
268  + F(i ,j-1,k+1,ncomp) - F(i ,j-1,k ,ncomp))
269  + gammaz * (F(i+1,j+1,k+1,ncomp) - F(i+1,j+1,k ,ncomp)
270  + F(i-1,j+1,k+1,ncomp) - F(i-1,j+1,k ,ncomp)
271  + F(i+1,j-1,k+1,ncomp) - F(i+1,j-1,k ,ncomp)
272  + F(i-1,j-1,k+1,ncomp) - F(i-1,j-1,k ,ncomp));
273 #elif (defined WARPX_DIM_XZ)
274  return alphaz * (F(i ,j+1,k ,ncomp) - F(i ,j ,k ,ncomp))
275  + betazx * (F(i+1,j+1,k ,ncomp) - F(i+1,j ,k ,ncomp)
276  + F(i-1,j+1,k ,ncomp) - F(i-1,j ,k ,ncomp));
277 #elif (defined WARPX_DIM_1D_Z)
278  return alphaz * (F(i+1 ,j,k ,ncomp) - F(i ,j ,k ,ncomp));
279 #endif
280  }
281 
284  template< typename T_Field>
286  static amrex::Real DownwardDz (
287  T_Field const& F,
288  amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
289  int const i, int const j, int const k, int const ncomp=0) {
290 
291  amrex::Real const inv_dz = coefs_z[0];
292 #if defined WARPX_DIM_3D
293  return inv_dz*( F(i,j,k,ncomp) - F(i,j,k-1,ncomp) );
294 #elif (defined WARPX_DIM_XZ)
295  return inv_dz*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) );
296 #elif (defined WARPX_DIM_1D_Z)
297  return inv_dz*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) );
298 #endif
299  }
300 
301 };
302 
303 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
#define AMREX_D_DECL(a, b, c)
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
i
Definition: check_interp_points_and_weights.py:174
cell_size
Definition: compute_domain.py:37
tuple dx
lab frame
Definition: stencil.py:429
beta
Definition: stencil.py:434
Definition: CartesianCKCAlgorithm.H:26
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real UpwardDz(amrex::Array4< amrex::Real > const &F, amrex::Real const *const coefs_z, int const n_coefs_z, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianCKCAlgorithm.H:246
static amrex::IntVect GetMaxGuardCell()
Returns maximum number of guard cells required by the field-solve.
Definition: CartesianCKCAlgorithm.H:123
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real UpwardDy(amrex::Array4< amrex::Real > const &F, amrex::Real const *const coefs_y, int const n_coefs_y, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianCKCAlgorithm.H:188
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real DownwardDx(T_Field const &F, amrex::Real const *const coefs_x, int const, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianCKCAlgorithm.H:170
static amrex::Real ComputeMaxDt(amrex::Real const *const dx)
Definition: CartesianCKCAlgorithm.H:107
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real UpwardDx(amrex::Array4< amrex::Real > const &F, amrex::Real const *const coefs_x, int const, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianCKCAlgorithm.H:131
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real DownwardDz(T_Field const &F, amrex::Real const *const coefs_z, int const, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianCKCAlgorithm.H:286
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real DownwardDy(T_Field const &F, amrex::Real const *const coefs_y, int const n_coefs_y, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianCKCAlgorithm.H:223
static void InitializeStencilCoefficients(std::array< amrex::Real, 3 > &cell_size, amrex::Vector< amrex::Real > &stencil_coefs_x, amrex::Vector< amrex::Real > &stencil_coefs_y, amrex::Vector< amrex::Real > &stencil_coefs_z)
Definition: CartesianCKCAlgorithm.H:28