WarpX
ShapeFactors.H
Go to the documentation of this file.
1 /* Copyright 2019-2021 Maxence Thevenet, Michael Rowan, Luca Fedeli, Axel Huebl
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef WARPX_SHAPEFACTORS_H_
8 #define WARPX_SHAPEFACTORS_H_
9 
10 #include "Utils/TextMsg.H"
11 
12 #include <AMReX.H>
13 #include <AMReX_GpuQualifiers.H>
14 
15 
27 template <int depos_order>
29 {
30  template< typename T >
33  T* const sx,
34  T xmid) const
35  {
36  if constexpr (depos_order == 0){
37  const auto j = static_cast<int>(xmid + T(0.5));
38  sx[0] = T(1.0);
39  return j;
40  }
41  else if constexpr (depos_order == 1){
42  const auto j = static_cast<int>(xmid);
43  const T xint = xmid - T(j);
44  sx[0] = T(1.0) - xint;
45  sx[1] = xint;
46  return j;
47  }
48  else if constexpr (depos_order == 2){
49  const auto j = static_cast<int>(xmid + T(0.5));
50  const T xint = xmid - T(j);
51  sx[0] = T(0.5)*(T(0.5) - xint)*(T(0.5) - xint);
52  sx[1] = T(0.75) - xint*xint;
53  sx[2] = T(0.5)*(T(0.5) + xint)*(T(0.5) + xint);
54  // index of the leftmost cell where particle deposits
55  return j-1;
56  }
57  else if constexpr (depos_order == 3){
58  const auto j = static_cast<int>(xmid);
59  const T xint = xmid - T(j);
60  sx[0] = (T(1.0))/(T(6.0))*(T(1.0) - xint)*(T(1.0) - xint)*(T(1.0) - xint);
61  sx[1] = (T(2.0))/(T(3.0)) - xint*xint*(T(1.0) - xint/(T(2.0)));
62  sx[2] = (T(2.0))/(T(3.0)) - (T(1.0) - xint)*(T(1.0) - xint)*(T(1.0) - T(0.5)*(T(1.0) - xint));
63  sx[3] = (T(1.0))/(T(6.0))*xint*xint*xint;
64  // index of the leftmost cell where particle deposits
65  return j-1;
66  }
67  else if constexpr (depos_order == 4){
68  const auto j = static_cast<int>(xmid + T(0.5));
69  const T xint = xmid - T(j);
70  sx[0] = (T(1.0))/(T(24.0))*(T(0.5) - xint)*(T(0.5) - xint)*(T(0.5) - xint)*(T(0.5) - xint);
71  sx[1] = (T(1.0))/(T(24.0))*(T(4.75) - T(11.0)*xint + T(4.0)*xint*xint*(T(1.5) + xint - xint*xint));
72  sx[2] = (T(1.0))/(T(24.0))*(T(14.375) + T(6.0)*xint*xint*(xint*xint - T(2.5)));
73  sx[3] = (T(1.0))/(T(24.0))*(T(4.75) + T(11.0)*xint + T(4.0)*xint*xint*(T(1.5) - xint - xint*xint));
74  sx[4] = (T(1.0))/(T(24.0))*(T(0.5) + xint)*(T(0.5) + xint)*(T(0.5) + xint)*(T(0.5)+xint);
75  // index of the leftmost cell where particle deposits
76  return j-2;
77  }
78  else{
79  WARPX_ABORT_WITH_MESSAGE("Unknown particle shape selected in Compute_shape_factor");
80  amrex::ignore_unused(sx, xmid);
81  }
82  return 0;
83  }
84 };
85 
86 
87 
93 template <int depos_order>
95 {
96  template< typename T >
99  T* const sx,
100  const T x_old,
101  const int i_new) const
102  {
103  if constexpr (depos_order == 0){
104  const auto i = static_cast<int>(std::floor(x_old + T(0.5)));
105  const int i_shift = i - i_new;
106  sx[1+i_shift] = T(1.0);
107  return i;
108  }
109  else if constexpr (depos_order == 1){
110  const auto i = static_cast<int>(std::floor(x_old));
111  const int i_shift = i - i_new;
112  const T xint = x_old - T(i);
113  sx[1+i_shift] = T(1.0) - xint;
114  sx[2+i_shift] = xint;
115  return i;
116  }
117  else if constexpr (depos_order == 2){
118  const auto i = static_cast<int>(x_old + T(0.5));
119  const int i_shift = i - (i_new + 1);
120  const T xint = x_old - T(i);
121  sx[1+i_shift] = T(0.5)*(T(0.5) - xint)*(T(0.5) - xint);
122  sx[2+i_shift] = T(0.75) - xint*xint;
123  sx[3+i_shift] = T(0.5)*(T(0.5) + xint)*(T(0.5) + xint);
124  // index of the leftmost cell where particle deposits
125  return i - 1;
126  }
127  else if constexpr (depos_order == 3){
128  const auto i = static_cast<int>(x_old);
129  const int i_shift = i - (i_new + 1);
130  const T xint = x_old - T(i);
131  sx[1+i_shift] = (T(1.0))/(T(6.0))*(T(1.0) - xint)*(T(1.0) - xint)*(T(1.0) - xint);
132  sx[2+i_shift] = (T(2.0))/(T(3.0)) - xint*xint*(T(1.0) - xint/(T(2.0)));
133  sx[3+i_shift] = (T(2.0))/(T(3.0)) - (T(1.0) - xint)*(T(1.0) - xint)*(T(1.0) - T(0.5)*(T(1.0) - xint));
134  sx[4+i_shift] = (T(1.0))/(T(6.0))*xint*xint*xint;
135  // index of the leftmost cell where particle deposits
136  return i - 1;
137  }
138  else if constexpr (depos_order == 4){
139  const auto i = static_cast<int>(x_old + T(0.5));
140  const int i_shift = i - (i_new + 2);
141  const T xint = x_old - T(i);
142  sx[1+i_shift] = (T(1.0))/(T(24.0))*(T(0.5) - xint)*(T(0.5) - xint)*(T(0.5) - xint)*(T(0.5) - xint);
143  sx[2+i_shift] = (T(1.0))/(T(24.0))*(T(4.75) - T(11.0)*xint + T(4.0)*xint*xint*(T(1.5) + xint - xint*xint));
144  sx[3+i_shift] = (T(1.0))/(T(24.0))*(T(14.375) + T(6.0)*xint*xint*(xint*xint - T(2.5)));
145  sx[4+i_shift] = (T(1.0))/(T(24.0))*(T(4.75) + T(11.0)*xint + T(4.0)*xint*xint*(T(1.5) - xint - xint*xint));
146  sx[5+i_shift] = (T(1.0))/(T(24.0))*(T(0.5) + xint)*(T(0.5) + xint)*(T(0.5) + xint)*(T(0.5)+xint);
147  // index of the leftmost cell where particle deposits
148  return i - 2;
149  }
150  else{
151  WARPX_ABORT_WITH_MESSAGE("Unknown particle shape selected in Compute_shifted_shape_factor");
152  amrex::ignore_unused(sx, x_old, i_new);
153  }
154  return 0;
155  }
156 };
157 
166 template <int depos_order>
168 {
169  template< typename T >
172  T* const sx_old,
173  T* const sx_new,
174  T xold,
175  T xnew) const
176  {
177  const T xmid = T(0.5)*(xnew + xold);
178  if constexpr (depos_order == 1){
179  const auto j = static_cast<int>(xmid);
180  const T xint_old = xold - T(j);
181  sx_old[0] = T(1.0) - xint_old;
182  sx_old[1] = xint_old;
183  //
184  const T xint_new = xnew - T(j);
185  sx_new[0] = T(1.0) - xint_new;
186  sx_new[1] = xint_new;
187  return j;
188  }
189  else if constexpr (depos_order == 2){
190  const auto j = static_cast<int>(xmid + T(0.5));
191  const T xint_old = xold - T(j);
192  sx_old[0] = T(0.5)*(T(0.5) - xint_old)*(T(0.5) - xint_old);
193  sx_old[1] = T(0.75) - xint_old*xint_old;
194  sx_old[2] = T(0.5)*(T(0.5) + xint_old)*(T(0.5) + xint_old);
195  //
196  const T xint_new = xnew - T(j);
197  sx_new[0] = T(0.5)*(T(0.5) - xint_new)*(T(0.5) - xint_new);
198  sx_new[1] = T(0.75) - xint_new*xint_new;
199  sx_new[2] = T(0.5)*(T(0.5) + xint_new)*(T(0.5) + xint_new);
200  // index of the leftmost cell where particle deposits
201  return j-1;
202  }
203  else if constexpr (depos_order == 3){
204  const auto j = static_cast<int>(xmid);
205  const T xint_old = xold - T(j);
206  sx_old[0] = T(1.0)/T(6.0)*(T(1.0) - xint_old)*(T(1.0) - xint_old)*(T(1.0) - xint_old);
207  sx_old[1] = T(2.0)/T(3.0) - xint_old*xint_old*(T(1.0) - xint_old/(T(2.0)));
208  sx_old[2] = T(2.0)/T(3.0) - (T(1.0) - xint_old)*(T(1.0) - xint_old)*(T(1.0) - T(0.5)*(T(1.0) - xint_old));
209  sx_old[3] = T(1.0)/T(6.0)*xint_old*xint_old*xint_old;
210  //
211  const T xint_new = xnew - T(j);
212  sx_new[0] = T(1.0)/T(6.0)*(T(1.0) - xint_new)*(T(1.0) - xint_new)*(T(1.0) - xint_new);
213  sx_new[1] = T(2.0)/T(3.0) - xint_new*xint_new*(T(1.0) - xint_new/(T(2.0)));
214  sx_new[2] = T(2.0)/T(3.0) - (T(1.0) - xint_new)*(T(1.0) - xint_new)*(T(1.0) - T(0.5)*(T(1.0) - xint_new));
215  sx_new[3] = T(1.0)/T(6.0)*xint_new*xint_new*xint_new;
216  // index of the leftmost cell where particle deposits
217  return j-1;
218  }
219  else if constexpr (depos_order == 4){
220  const auto j = static_cast<int>(xmid + T(0.5));
221  const T xint_old = xold - T(j);
222  sx_old[0] = (T(1.0))/(T(24.0))*(T(0.5) - xint_old)*(T(0.5) - xint_old)*(T(0.5) - xint_old)*(T(0.5) - xint_old);
223  sx_old[1] = (T(1.0))/(T(24.0))*(T(4.75) - T(11.0)*xint_old + T(4.0)*xint_old*xint_old*(T(1.5) + xint_old - xint_old*xint_old));
224  sx_old[2] = (T(1.0))/(T(24.0))*(T(14.375) + T(6.0)*xint_old*xint_old*(xint_old*xint_old - T(2.5)));
225  sx_old[3] = (T(1.0))/(T(24.0))*(T(4.75) + T(11.0)*xint_old + T(4.0)*xint_old*xint_old*(T(1.5) - xint_old - xint_old*xint_old));
226  sx_old[4] = (T(1.0))/(T(24.0))*(T(0.5) + xint_old)*(T(0.5) + xint_old)*(T(0.5) + xint_old)*(T(0.5)+xint_old);
227  //
228  const T xint_new = xnew - T(j);
229  sx_new[0] = (T(1.0))/(T(24.0))*(T(0.5) - xint_new)*(T(0.5) - xint_new)*(T(0.5) - xint_new)*(T(0.5) - xint_new);
230  sx_new[1] = (T(1.0))/(T(24.0))*(T(4.75) - T(11.0)*xint_new + T(4.0)*xint_new*xint_new*(T(1.5) + xint_new - xint_new*xint_new));
231  sx_new[2] = (T(1.0))/(T(24.0))*(T(14.375) + T(6.0)*xint_new*xint_new*(xint_new*xint_new - T(2.5)));
232  sx_new[3] = (T(1.0))/(T(24.0))*(T(4.75) + T(11.0)*xint_new + T(4.0)*xint_new*xint_new*(T(1.5) - xint_new - xint_new*xint_new));
233  sx_new[4] = (T(1.0))/(T(24.0))*(T(0.5) + xint_new)*(T(0.5) + xint_new)*(T(0.5) + xint_new)*(T(0.5)+xint_new);
234  //
235  // index of the leftmost cell where particle deposits
236  return j-2;
237  }
238  else{
239  WARPX_ABORT_WITH_MESSAGE("Unknown particle shape selected in Compute_shape_factor_pair");
240  amrex::ignore_unused(sx_old, sx_new, xold, xnew);
241  }
242  return 0;
243  }
244 };
245 
246 #endif // WARPX_SHAPEFACTORS_H_
#define AMREX_FORCE_INLINE
#define AMREX_GPU_HOST_DEVICE
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition: TextMsg.H:15
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
i
Definition: check_interp_points_and_weights.py:174
Definition: ShapeFactors.H:168
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int operator()(T *const sx_old, T *const sx_new, T xold, T xnew) const
Definition: ShapeFactors.H:171
Definition: ShapeFactors.H:29
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int operator()(T *const sx, T xmid) const
Definition: ShapeFactors.H:32
Definition: ShapeFactors.H:95
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int operator()(T *const sx, const T x_old, const int i_new) const
Definition: ShapeFactors.H:98