WarpX
CylindricalYeeAlgorithm.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_CYLINDRICAL_YEE_H_
9 #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_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 <array>
19 #include <cmath>
20 
21 
27 
29  std::array<amrex::Real,3>& cell_size,
30  amrex::Vector<amrex::Real>& stencil_coefs_r,
31  amrex::Vector<amrex::Real>& stencil_coefs_z ) {
32 
33  using namespace amrex;
34  // Store the inverse cell size along each direction in the coefficients
35  stencil_coefs_r.resize(1);
36  stencil_coefs_r[0] = 1._rt/cell_size[0]; // 1./dr
37  stencil_coefs_z.resize(1);
38  stencil_coefs_z[0] = 1._rt/cell_size[2]; // 1./dz
39  }
40 
46  static amrex::Real ComputeMaxDt ( amrex::Real const * const dx,
47  int const n_rz_azimuthal_modes ) {
48  using namespace amrex::literals;
49  // In the rz case, the Courant limit has been evaluated
50  // semi-analytically by R. Lehe, and resulted in the following
51  // coefficients.
52  std::array< amrex::Real, 6 > const multimode_coeffs = {{ 0.2105_rt, 1.0_rt, 3.5234_rt, 8.5104_rt, 15.5059_rt, 24.5037_rt }};
53  const amrex::Real multimode_alpha = (n_rz_azimuthal_modes < 7)?
54  multimode_coeffs[n_rz_azimuthal_modes-1]: // Use the table of the coefficients
55  (n_rz_azimuthal_modes - 1._rt)*(n_rz_azimuthal_modes - 1._rt) - 0.4_rt; // Use a realistic extrapolation
56  const amrex::Real delta_t = 1._rt / ( std::sqrt(
57  (1._rt + multimode_alpha) / (dx[0]*dx[0])
58  + 1._rt / (dx[1]*dx[1])
59  ) * PhysConst::c );
60  return delta_t;
61  }
62 
67  // The cylindrical solver requires one guard cell in each dimension
68  return amrex::IntVect{AMREX_D_DECL(1,1,1)};
69  }
70 
76  static amrex::Real UpwardDrr_over_r (
78  amrex::Real const r, amrex::Real const dr,
79  amrex::Real const * const coefs_r, int const n_coefs_r,
80  int const i, int const j, int const k, int const comp ) {
81 
82  using namespace amrex;
83  ignore_unused(n_coefs_r);
84 
85  Real const inv_dr = coefs_r[0];
86  return 1._rt/r * inv_dr*( (r+0.5_rt*dr)*F(i+1,j,k,comp) - (r-0.5_rt*dr)*F(i,j,k,comp) );
87  }
88 
94  static amrex::Real DownwardDrr_over_r (
96  amrex::Real const r, amrex::Real const dr,
97  amrex::Real const * const coefs_r, int const n_coefs_r,
98  int const i, int const j, int const k, int const comp ) {
99 
100  using namespace amrex;
101  ignore_unused(n_coefs_r);
102 
103  Real const inv_dr = coefs_r[0];
104  return 1._rt/r * inv_dr*( (r+0.5_rt*dr)*F(i,j,k,comp) - (r-0.5_rt*dr)*F(i-1,j,k,comp) );
105  }
106 
110  static amrex::Real UpwardDr (
112  amrex::Real const * const coefs_r, int const n_coefs_r,
113  int const i, int const j, int const k, int const comp ) {
114 
115  using namespace amrex;
116  ignore_unused(n_coefs_r);
117 
118  Real const inv_dr = coefs_r[0];
119  return inv_dr*( F(i+1,j,k,comp) - F(i,j,k,comp) );
120  }
121 
125  static amrex::Real DownwardDr (
127  amrex::Real const * const coefs_r, int const n_coefs_r,
128  int const i, int const j, int const k, int const comp ) {
129 
130  using namespace amrex;
131  ignore_unused(n_coefs_r);
132 
133  Real const inv_dr = coefs_r[0];
134  return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) );
135  }
136 
140  template< typename T_Field>
142  static amrex::Real Dr_rDr_over_r (
143  T_Field const& F,
144  amrex::Real const r, amrex::Real const dr,
145  amrex::Real const * const coefs_r, int const /*n_coefs_r*/,
146  int const i, int const j, int const k, int const comp ) {
147 
148  using namespace amrex;
149 
150  Real const inv_dr2 = coefs_r[0]*coefs_r[0];
151  return 1._rt/r * inv_dr2*( (r+0.5_rt*dr)*(F(i+1,j,k,comp) - F(i,j,k,comp))
152  - (r-0.5_rt*dr)*(F(i,j,k,comp) - F(i-1,j,k,comp)) );
153  }
154 
158  static amrex::Real UpwardDz (
160  amrex::Real const * const coefs_z, int const n_coefs_z,
161  int const i, int const j, int const k, int const comp ) {
162 
163  amrex::ignore_unused(n_coefs_z);
164 
165  amrex::Real const inv_dz = coefs_z[0];
166  return inv_dz*( F(i,j+1,k,comp) - F(i,j,k,comp) );
167  }
168 
172  static amrex::Real DownwardDz (
174  amrex::Real const * const coefs_z, int const n_coefs_z,
175  int const i, int const j, int const k, int const comp ) {
176 
177  amrex::ignore_unused(n_coefs_z);
178 
179  amrex::Real const inv_dz = coefs_z[0];
180  return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) );
181  }
182 
185  template< typename T_Field>
187  static amrex::Real Dzz (
188  T_Field const& F,
189  amrex::Real const * const coefs_z, int const /*n_coefs_z*/,
190  int const i, int const j, int const k, int const ncomp=0 ) {
191 
192  using namespace amrex;
193  Real const inv_dz2 = coefs_z[0]*coefs_z[0];
194 
195  return inv_dz2*( F(i,j-1,k,ncomp) - 2._rt*F(i,j,k,ncomp) + F(i,j+1,k,ncomp) );
196  }
197 
198 };
199 
200 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CYLINDRICAL_YEE_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
Definition: CylindricalYeeAlgorithm.H:26
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: CylindricalYeeAlgorithm.H:187
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 comp)
Definition: CylindricalYeeAlgorithm.H:158
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real UpwardDr(amrex::Array4< amrex::Real > const &F, amrex::Real const *const coefs_r, int const n_coefs_r, int const i, int const j, int const k, int const comp)
Definition: CylindricalYeeAlgorithm.H:110
static amrex::IntVect GetMaxGuardCell()
Returns maximum number of guard cells required by the field-solve.
Definition: CylindricalYeeAlgorithm.H:66
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real Dr_rDr_over_r(T_Field const &F, amrex::Real const r, amrex::Real const dr, amrex::Real const *const coefs_r, int const, int const i, int const j, int const k, int const comp)
Definition: CylindricalYeeAlgorithm.H:142
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real DownwardDr(amrex::Array4< amrex::Real > const &F, amrex::Real const *const coefs_r, int const n_coefs_r, int const i, int const j, int const k, int const comp)
Definition: CylindricalYeeAlgorithm.H:125
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real DownwardDrr_over_r(amrex::Array4< amrex::Real > const &F, amrex::Real const r, amrex::Real const dr, amrex::Real const *const coefs_r, int const n_coefs_r, int const i, int const j, int const k, int const comp)
Definition: CylindricalYeeAlgorithm.H:94
static amrex::Real ComputeMaxDt(amrex::Real const *const dx, int const n_rz_azimuthal_modes)
Definition: CylindricalYeeAlgorithm.H:46
static void InitializeStencilCoefficients(std::array< amrex::Real, 3 > &cell_size, amrex::Vector< amrex::Real > &stencil_coefs_r, amrex::Vector< amrex::Real > &stencil_coefs_z)
Definition: CylindricalYeeAlgorithm.H:28
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real DownwardDz(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 comp)
Definition: CylindricalYeeAlgorithm.H:172
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real UpwardDrr_over_r(amrex::Array4< amrex::Real > const &F, amrex::Real const r, amrex::Real const dr, amrex::Real const *const coefs_r, int const n_coefs_r, int const i, int const j, int const k, int const comp)
Definition: CylindricalYeeAlgorithm.H:76