WarpX
CartesianYeeAlgorithm.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_YEE_H_
9 #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_
10 
11 #include "FieldAccessorFunctors.H"
12 #include "Utils/WarpXConst.H"
13 
14 #include <AMReX.H>
15 #include <AMReX_Array4.H>
16 #include <AMReX_Gpu.H>
17 #include <AMReX_REAL.H>
18 
19 #include <array>
20 #include <cmath>
21 
22 
28 
30  std::array<amrex::Real,3>& cell_size,
31  amrex::Vector<amrex::Real>& stencil_coefs_x,
32  amrex::Vector<amrex::Real>& stencil_coefs_y,
33  amrex::Vector<amrex::Real>& stencil_coefs_z ) {
34 
35  using namespace amrex;
36  // Store the inverse cell size along each direction in the coefficients
37  stencil_coefs_x.resize(1);
38  stencil_coefs_x[0] = 1._rt/cell_size[0];
39  stencil_coefs_y.resize(1);
40  stencil_coefs_y[0] = 1._rt/cell_size[1];
41  stencil_coefs_z.resize(1);
42  stencil_coefs_z[0] = 1._rt/cell_size[2];
43  }
44 
48  static amrex::Real ComputeMaxDt ( amrex::Real const * const dx ) {
49  using namespace amrex::literals;
50  amrex::Real const delta_t = 1._rt / ( std::sqrt( AMREX_D_TERM(
51  1._rt / (dx[0]*dx[0]),
52  + 1._rt / (dx[1]*dx[1]),
53  + 1._rt / (dx[2]*dx[2])
54  ) ) * PhysConst::c );
55  return delta_t;
56  }
57 
62  // The yee solver requires one guard cell in each dimension
63  return amrex::IntVect{AMREX_D_DECL(1,1,1)};
64  }
65 
69  static amrex::Real UpwardDx (
71  amrex::Real const * const coefs_x, int const /*n_coefs_x*/,
72  int const i, int const j, int const k, int const ncomp=0 ) {
73 
74  using namespace amrex;
75 #if (defined WARPX_DIM_1D_Z)
76  amrex::ignore_unused(F, coefs_x, i, j, k, ncomp);
77  return 0._rt; // 1D Cartesian: derivative along x is 0
78 #else
79  amrex::Real const inv_dx = coefs_x[0];
80  return inv_dx*( F(i+1,j,k,ncomp) - F(i,j,k,ncomp) );
81 #endif
82  }
83 
86  template< typename T_Field>
88  static amrex::Real DownwardDx (
89  T_Field const& F,
90  amrex::Real const * const coefs_x, int const /*n_coefs_x*/,
91  int const i, int const j, int const k, int const ncomp=0 ) {
92 
93  using namespace amrex;
94 #if (defined WARPX_DIM_1D_Z)
95  amrex::ignore_unused(F, coefs_x, i, j, k, ncomp);
96  return 0._rt; // 1D Cartesian: derivative along x is 0
97 #else
98  amrex::Real const inv_dx = coefs_x[0];
99  return inv_dx*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) );
100 #endif
101  }
102 
105  template< typename T_Field>
107  static amrex::Real Dxx (
108  T_Field const& F,
109  amrex::Real const * const coefs_x, int const /*n_coefs_x*/,
110  int const i, int const j, int const k, int const ncomp=0 ) {
111 
112  using namespace amrex;
113 #if (defined WARPX_DIM_1D_Z)
114  amrex::ignore_unused(F, coefs_x, i, j, k, ncomp);
115  return 0._rt; // 1D Cartesian: derivative along x is 0
116 #else
117  amrex::Real const inv_dx2 = coefs_x[0]*coefs_x[0];
118  return inv_dx2*( F(i-1,j,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i+1,j,k,ncomp) );
119 #endif
120  }
121 
125  static amrex::Real UpwardDy (
127  amrex::Real const * const coefs_y, int const n_coefs_y,
128  int const i, int const j, int const k, int const ncomp=0 ) {
129 
130  using namespace amrex;
131 #if defined WARPX_DIM_3D
132  Real const inv_dy = coefs_y[0];
133  amrex::ignore_unused(n_coefs_y);
134  return inv_dy*( F(i,j+1,k,ncomp) - F(i,j,k,ncomp) );
135 #elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z)
136  amrex::ignore_unused(F, coefs_y, n_coefs_y,
137  i, j, k, ncomp);
138  return 0._rt; // 1D and 2D Cartesian: derivative along y is 0
139 #else
140  amrex::ignore_unused(F, coefs_y, n_coefs_y,
141  i, j, k, ncomp);
142 #endif
143  }
144 
147  template< typename T_Field>
149  static amrex::Real DownwardDy (
150  T_Field const& F,
151  amrex::Real const * const coefs_y, int const n_coefs_y,
152  int const i, int const j, int const k, int const ncomp=0 ) {
153 
154  using namespace amrex;
155 #if defined WARPX_DIM_3D
156  Real const inv_dy = coefs_y[0];
157  amrex::ignore_unused(n_coefs_y);
158  return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) );
159 #elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z)
160  amrex::ignore_unused(F, coefs_y, n_coefs_y,
161  i, j, k, ncomp);
162  return 0._rt; // 1D and 2D Cartesian: derivative along y is 0
163 #else
164  amrex::ignore_unused(F, coefs_y, n_coefs_y,
165  i, j, k, ncomp);
166 #endif
167  }
168 
171  template< typename T_Field>
173  static amrex::Real Dyy (
174  T_Field const& F,
175  amrex::Real const * const coefs_y, int const /*n_coefs_y*/,
176  int const i, int const j, int const k, int const ncomp=0 ) {
177 
178  using namespace amrex;
179 #if defined WARPX_DIM_3D
180  Real const inv_dy2 = coefs_y[0]*coefs_y[0];
181  return inv_dy2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) );
182 #elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z)
183  amrex::ignore_unused(F, coefs_y, i, j, k, ncomp);
184  return 0._rt; // 1D and 2D Cartesian: derivative along y is 0
185 #endif
186  }
187 
191  static amrex::Real UpwardDz (
193  amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
194  int const i, int const j, int const k, int const ncomp=0 ) {
195 
196  using namespace amrex;
197  Real const inv_dz = coefs_z[0];
198 #if defined WARPX_DIM_3D
199  return inv_dz*( F(i,j,k+1,ncomp) - F(i,j,k,ncomp) );
200 #elif (defined WARPX_DIM_XZ)
201  return inv_dz*( F(i,j+1,k,ncomp) - F(i,j,k,ncomp) );
202 #elif (defined WARPX_DIM_1D_Z)
203  return inv_dz*( F(i+1,j,k,ncomp) - F(i,j,k,ncomp) );
204 #endif
205  }
206 
209  template< typename T_Field>
211  static amrex::Real DownwardDz (
212  T_Field const& F,
213  amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
214  int const i, int const j, int const k, int const ncomp=0 ) {
215 
216  using namespace amrex;
217  Real const inv_dz = coefs_z[0];
218 #if defined WARPX_DIM_3D
219  return inv_dz*( F(i,j,k,ncomp) - F(i,j,k-1,ncomp) );
220 #elif (defined WARPX_DIM_XZ)
221  return inv_dz*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) );
222 #elif (defined WARPX_DIM_1D_Z)
223  return inv_dz*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) );
224 #endif
225  }
226 
229  template< typename T_Field>
231  static amrex::Real Dzz (
232  T_Field const& F,
233  amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
234  int const i, int const j, int const k, int const ncomp=0 ) {
235 
236  using namespace amrex;
237  Real const inv_dz2 = coefs_z[0]*coefs_z[0];
238 #if defined WARPX_DIM_3D
239  return inv_dz2*( F(i,j,k-1,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j,k+1,ncomp) );
240 #elif (defined WARPX_DIM_XZ)
241  return inv_dz2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) );
242 #elif (defined WARPX_DIM_1D_Z)
243  return inv_dz2*( F(i-1,j,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i+1,j,k,ncomp) );
244 #endif
245  }
246 
247 };
248 
249 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
#define AMREX_D_TERM(a, b, c)
#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
Definition: CartesianYeeAlgorithm.H:27
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real Dxx(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: CartesianYeeAlgorithm.H:107
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: CartesianYeeAlgorithm.H:88
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, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianYeeAlgorithm.H:191
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: CartesianYeeAlgorithm.H:125
static amrex::IntVect GetMaxGuardCell()
Returns maximum number of guard cells required by the field-solve.
Definition: CartesianYeeAlgorithm.H:61
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: CartesianYeeAlgorithm.H:69
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: CartesianYeeAlgorithm.H:149
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real Dzz(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: CartesianYeeAlgorithm.H:231
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: CartesianYeeAlgorithm.H:211
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real Dyy(T_Field const &F, amrex::Real const *const coefs_y, int const, int const i, int const j, int const k, int const ncomp=0)
Definition: CartesianYeeAlgorithm.H:173
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: CartesianYeeAlgorithm.H:29
static amrex::Real ComputeMaxDt(amrex::Real const *const dx)
Definition: CartesianYeeAlgorithm.H:48