8 #ifndef WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_
9 #define WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_CKC_H_
34 using namespace amrex;
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;
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;
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;
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;
108 #if (defined WARPX_DIM_1D_Z)
110 #elif (defined WARPX_DIM_XZ)
115 amrex::Real
const delta_t = std::min(
dx[0], std::min(
dx[1],
dx[2] ) ) /
PhysConst::c;
133 amrex::Real
const *
const coefs_x,
int const ,
134 int const i,
int const j,
int const k,
int const ncomp=0 ) {
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];
141 #if defined WARPX_DIM_3D
142 amrex::Real
const betaxy = coefs_x[2];
143 amrex::Real
const gammax = coefs_x[4];
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)
168 template<
typename T_Field>
172 amrex::Real
const *
const coefs_x,
int const ,
173 int const i,
int const j,
int const k,
int const ncomp=0 ) {
175 using namespace amrex;
176 #if (defined WARPX_DIM_1D_Z)
180 amrex::Real
const inv_dx = coefs_x[0];
181 return inv_dx*( F(
i,j,k,ncomp) - F(
i-1,j,k,ncomp) );
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 ) {
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];
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)
221 template<
typename T_Field>
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 ) {
228 using namespace amrex;
229 #if defined WARPX_DIM_3D
230 amrex::Real
const inv_dy = coefs_y[0];
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)
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 ) {
251 using namespace amrex;
255 Real
const alphaz = coefs_z[1];
256 #if (defined WARPX_DIM_3D || WARPX_DIM_XZ)
257 Real
const betazx = coefs_z[2];
259 #if defined WARPX_DIM_3D
260 Real
const betazy = coefs_z[3];
261 Real
const gammaz = coefs_z[4];
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));
284 template<
typename T_Field>
288 amrex::Real
const *
const coefs_z,
int const ,
289 int const i,
int const j,
int const k,
int const ncomp=0) {
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) );
#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