XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
get_potential_grad.hpp
Go to the documentation of this file.
1 #ifndef GET_POTENTIAL_GRAD
2 #define GET_POTENTIAL_GRAD
3 
4 #include "sml.hpp"
5 #include "grid.hpp"
7 #include "magnetic_field.hpp"
8 #include "my_mirror_view.hpp"
9 #include "gradparx2.hpp"
11 #include "gradient_matrices.hpp"
12 #include "gyro_avg_mat.hpp"
13 #include "task_group.hpp"
14 #include "grid_field.hpp"
15 #include "grid_deriv.hpp"
16 
17 // Fortran pointer retrievals
18 
19 extern "C" rtype* get_psn_pbd0_2_iseg_loc();
20 extern "C" int get_psn_pbd0_2_nseg();
21 
22 // Fortran routines
23 
24 extern "C" void get_pot_epar_em_filter(double* tmp, double* E_para_em);
25 
26 
27 template<class Device>
29  View<double**,CLayout,Device> field;
30 
31  public:
32 
33  AlternatingStorage(const std::string& name, int nnode)
34  : field(NoInit(name), 2, nnode)
35  {}
36 
37  // Functions to access the "alternating" Views
38  View<double*,CLayout, Device> left(int i_plane){
39  int i_left_alt = i_plane%2;
40  return my_subview(field,i_left_alt);
41  }
42 
43  View<double*,CLayout, Device> right(int i_plane){
44  int i_right_alt = (i_plane+1)%2;
45  return my_subview(field,i_right_alt);
46  }
47 };
48 
49 template<class Device, class DeviceOut>
51 
52  int nnode;
53  int nrho;
54  int nphi;
55  int ndim;
56 
57  View<double**,CLayout,Device> input_potential;
61  View<double**,CLayout,Device> potential;
62  View<double***,CLayout, Device> gradient;
63  View<double****,CLayout, Device> gradient_rho;
64  View<double***,CLayout, Device> potential_rho;
65 
66  View<double*,CLayout,Device> scratch;
67 
68  // Additional objects that should go elsewhere eventually, but are here
69  // because they are needed for get_pot_grad and can be reused for the different fields
74 
75  GetPotentialGradTemp(const Simulation<DeviceType>& sml, const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const MagneticField<DeviceType>& magnetic_field, int n_input_potential_planes, bool gyroavg_requested, const GyroAverageMatrices<HostType>& gyro_avg_matrices_in = GyroAverageMatrices<HostType>())
76  : nnode(grid.nnode),
77  nrho(gyroavg_requested ? gyro_avg_matrices_in.nrho : 0),
78  nphi((PIT_GLOBAL==PhiInterpType::Planes) ? 2 : 1),
80  potential_alt("potential_alt", nnode),
81  gradient_r_alt("gradient_r_alt", nnode),
82  gradient_z_alt("gradient_r_alt", nnode),
83  potential(NoInit("potential"), nphi, nnode),
84  gradient(NoInit("gradient"), ndim, nphi, nnode),
85  potential_rho(NoInit("potential_rho"), nrho+1, nphi, nnode),
86  gradient_rho(NoInit("gradient_rho"), nrho+1, nphi, nnode, ndim),
87  gyro_avg_matrices(gyro_avg_matrices_in.template mirror<DeviceType>()) // Send gyromatrices to GPU
88  {
89  // Use scratch for transpose on device if the fields can't be written directly
90  if(!std::is_same<DeviceOut,Device>()){
91  scratch = View<double*,CLayout,Device>(NoInit("scratch"), gradient_rho.size());
92  input_potential = View<double**,CLayout,Device>(NoInit("input_potential"), n_input_potential_planes, nnode);
93  }
94 
95  GPTLstart("GET_POT_GRAD_MAT_SETUP");
96 #ifdef NO_FORTRAN_MODULES
97  grad_matrices = grid.gradient_matrices_h.template mirror<DeviceType>();
98 #else
99  grad_matrices = GradientMatrices<DeviceType>(true); // Copy in fortran data
100 #endif
101  GPTLstop("GET_POT_GRAD_MAT_SETUP");
102 
103  // This could be set up only once (if it doesnt take up too much device memory)
104  GPTLstart("GET_POT_GRAD_FF_SETUP");
106  GPTLstop("GET_POT_GRAD_FF_SETUP");
107 
108  // Set up for parallel gradient
109  GPTLstart("GET_POT_GRAD_GPTX_SETUP");
110 #ifdef NO_FORTRAN_MODULES
111  if(!sml.is_XGCa) gptx = GradParXTmp(grid, magnetic_field.bt_sign);
112 #else
113  if(!sml.is_XGCa) gptx = GradParXTmp(grid, get_psn_pbd0_2_nseg(), magnetic_field.bt_sign, get_psn_pbd0_2_iseg_loc());
114 #endif
115  GPTLstop("GET_POT_GRAD_GPTX_SETUP");
116  }
117 };
118 
119 template<class Device>
120 struct Field00{
123 
124  View<Field<VarType::Scalar,PhiInterpType::None>*,CLayout, Device> pot_managed; // Use if device != host
125  View<double*,CLayout, Device, Kokkos::MemoryTraits<Kokkos::Unmanaged>> pot;
126  View<double*,CLayout, Device> r;
127  View<double*,CLayout, Device> z;
128 
130  : is_provided(false)
131  {}
132 
133  template<class DeviceIn>
135  : r(NoInit("r"), field_in.f.extent(0)),
136  z(NoInit("z"), field_in.f.extent(0)),
137  is_provided(true),
138  discard_when_basis_is_one(discard_when_basis_is_one_in)
139  {
140  // Create a mirror if Device is different from the input DeviceIn
141  pot_managed = my_mirror_view(field_in.f, Device());
142  mirror_copy(pot_managed, field_in.f);
143  pot = View<double*,CLayout,Device, Kokkos::MemoryTraits<Kokkos::Unmanaged>>((double*)(pot_managed.data()), pot_managed.layout());
144  }
145 
146  void calculate_gradient(const Grid<DeviceType>& grid, const GradientMatrices<DeviceType>& grad_matrices){
147  grid_deriv(grad_matrices, grid,pot,r,z, discard_when_basis_is_one);
148 
149  auto r_tmp = r; // Use these to avoid using member function in kernel
150  auto z_tmp = z; // Use these to avoid using member function in kernel
151 
152  // rh This is general now --->
153  // We compute the electric field --> flip sign
154  // E=-E
155  Kokkos::parallel_for("get_pot_grad_flip_00", Kokkos::RangePolicy<ExSpace>(0,grid.nnode), KOKKOS_LAMBDA( const int i ){
156  r_tmp(i)=-r_tmp(i);
157  z_tmp(i)=-z_tmp(i);
158  });
159  }
160 
161  template<class DeviceOut, PhiInterpType PIT>
162  void set_field(const Grid<DeviceType>& grid, const FieldFollowingCoordinates& ff, View<Field<VarType::Vector2D,PIT>*,CLayout,DeviceOut>& field00_ff_h){
163  if (is_provided){
164  // Copy to two planes if field following is needed
165  int nphi = (PIT==PhiInterpType::Planes) ? 2 : 1;
166  View<double**,CLayout, Device> r_tmp(NoInit("r_tmp"), nphi, r.size());
167  View<double**,CLayout, Device> z_tmp(NoInit("z_tmp"), nphi, z.size());
168 
169  for(int iphi=0; iphi<nphi; iphi++){
170  Kokkos::deep_copy(my_subview(r_tmp, iphi), r);
171  Kokkos::deep_copy(my_subview(z_tmp, iphi), z);
172  }
173 
174  if(PIT==PhiInterpType::Planes){
175  ff.cnvt_grid_real2ff(grid,r_tmp);
176  ff.cnvt_grid_real2ff(grid,z_tmp);
177  }
178 
179  auto field00_ff = my_mirror_view(field00_ff_h, Device());
180  for(int iphi=0; iphi<nphi; iphi++){
181  Kokkos::parallel_for("set_field", Kokkos::RangePolicy<ExSpace>(0,field00_ff_h.extent(0)), KOKKOS_LAMBDA( const int i ){
182 #ifdef XGC1
183  field00_ff(i).V[iphi][0]=r_tmp(iphi,i);
184  field00_ff(i).V[iphi][1]=z_tmp(iphi,i);
185 #else
186  field00_ff(i).E[0]=r_tmp(iphi,i);
187  field00_ff(i).E[1]=z_tmp(iphi,i);
188 #endif
189  });
190  }
191  Kokkos::fence();
192  // Copy if not a mirror
193  mirror_copy(field00_ff_h, field00_ff);
194  }else{
195  Kokkos::parallel_for("set_field", Kokkos::RangePolicy<typename DeviceOut::execution_space>(0,field00_ff_h.extent(0)), KOKKOS_LAMBDA( const int i ){
196 #ifdef XGC1
197  field00_ff_h(i).V[0][0]=0.0;
198  field00_ff_h(i).V[0][1]=0.0;
199  field00_ff_h(i).V[1][0]=0.0;
200  field00_ff_h(i).V[1][1]=0.0;
201 #else
202  field00_ff_h(i).E[0]=0.0;
203  field00_ff_h(i).E[1]=0.0;
204 #endif
205  });
206  }
207  }
208 };
209 
210 template<class Device, class DeviceOut>
211 struct EMParField{
212  bool requested;
213 
214  View<double**,CLayout, Device> field;
215  View<double***,CLayout, Device> field_rho;
216 
218 
220  : requested(false)
221  {}
222 
224  requested = true;
225  field_out = output_field;
226  field_rho = View<double***,CLayout, Device>(NoInit("E_para_em_rho"), field_out.f.extent(1), 2, field_out.f.extent(0));
227  field = View<double**,CLayout, Device>(NoInit("E_para_em"), 2, field_out.f.extent(0));
228  }
229 
230  void calculate(const Grid<DeviceType>& grid, const View<double**,CLayout, Device>& field_para){
231  // Obtain E_para with the same filters that are applied in push_As
232  // This filtered E_para must be used in the equation of motion
233  // for dA_s/dt in case of pullback-mode 4.
234  GPTLstart("GET_POT_EPAR_EM");
235  // tmp_copy gets copied in and modified internally
236  View<double**,CLayout, HostType> tmp_copy(NoInit("tmp"),2,grid.nnode);
237  Kokkos::deep_copy(tmp_copy, field_para);
238 
239  // Output:
240  auto field_h = my_mirror_view(field, HostType());
241 
242  get_pot_epar_em_filter(tmp_copy.data(), field_h.data());
243 
244  // Copy result of get_pot_epar_em_filter back to device
245  mirror_copy(field, field_h);
246  GPTLstop("GET_POT_EPAR_EM");
247  }
248 };
249 
250 template<class DeviceIn, class DeviceOut, VarType VT, PhiInterpType PIT, TorType TT, KinType KT>
252  View<double**,CLayout,DeviceIn, Kokkos::MemoryTraits<Kokkos::Unmanaged>> input_potential;
260 
261  // Constructor if passed an unmanaged view
262  GetPotGradFieldArgs(const View<double**,CLayout,DeviceIn, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& input_potential_in, bool ignore_poloidal_dpot_in=false)
263  : input_potential(input_potential_in),
264  ignore_poloidal_dpot(ignore_poloidal_dpot_in),
265  potential_is_requested(false),
266  gradient_is_requested(false)
267  {}
268 
269  // Constructor if passed a 1D managed view
270  GetPotGradFieldArgs(const View<double*,CLayout,DeviceIn>& input_potential_in, bool ignore_poloidal_dpot_in=false)
271  : input_potential(View<double**,CLayout,DeviceIn,Kokkos::MemoryTraits<Kokkos::Unmanaged>>((double*)(input_potential_in.data()), 1, input_potential_in.extent(0))),
272  ignore_poloidal_dpot(ignore_poloidal_dpot_in),
273  potential_is_requested(false),
274  gradient_is_requested(false)
275  {}
276 
278  potential = potential_in;
279  potential_is_requested = true;
280  }
281 
283  gradient = gradient_in;
284  gradient_is_requested = true;
285  }
286 
288  E_para_em.request(output_field);
289  }
290 };
291 
293 
294 template<class DeviceIn, class DeviceOut, VarType VT, PhiInterpType PIT, TorType TT, KinType KT>
295 void get_field_grad(const Grid<DeviceType>& grid,
296  const DomainDecomposition<DeviceType>& pol_decomp,
300 
301 #endif
VarType
Definition: field.hpp:16
AlternatingStorage(const std::string &name, int nnode)
Definition: get_potential_grad.hpp:33
View< double **, CLayout, Device > input_potential
Definition: get_potential_grad.hpp:57
Definition: field_following_coordinates.hpp:9
static int GPTLstart(const char *name)
Definition: timer_macro.hpp:9
int nrho
Definition: get_potential_grad.hpp:53
View< double *, CLayout, Device > left(int i_plane)
Definition: get_potential_grad.hpp:38
Definition: boundary.hpp:6
constexpr VarType vec2d_if_axisym()
View< double *, CLayout, Device > scratch
Definition: get_potential_grad.hpp:66
int ndim
Definition: get_potential_grad.hpp:55
View< Field< VarType::Scalar, PhiInterpType::None > *, CLayout, Device > pot_managed
Definition: get_potential_grad.hpp:124
void mirror_copy(T1 &view_dest, const T2 &view_src)
Definition: my_mirror_view.hpp:122
bool ignore_poloidal_dpot
Definition: get_potential_grad.hpp:255
void get_pot_epar_em_filter(double *tmp, double *E_para_em)
bool is_provided
Definition: get_potential_grad.hpp:121
Field00()
Definition: get_potential_grad.hpp:129
Definition: sml.hpp:8
static constexpr bool is_XGCa
Equivalent to the preprocessor flag for now.
Definition: sml.hpp:19
Kokkos::Device< HostExSpace, HostMemSpace > HostType
Definition: space_settings.hpp:56
FieldFollowingCoordinates ff
Definition: get_potential_grad.hpp:72
View< double **, CLayout, Device > potential
Definition: get_potential_grad.hpp:61
int nphi
Definition: get_potential_grad.hpp:54
Definition: magnetic_field.hpp:12
Definition: get_potential_grad.hpp:251
int get_psn_pbd0_2_nseg()
void request_para_em(const GridField< DeviceOut, VarType::Scalar, PhiInterpType::Planes, TorType::OnePlane, KinType::GyroKin > &output_field)
Definition: get_potential_grad.hpp:287
void cnvt_grid_real2ff(const Grid< DeviceType > &grid, const View< double **, CLayout, DeviceType > &input, const View< double **, Kokkos::LayoutRight, DeviceType > &output) const
Definition: field_following_coordinates.hpp:121
bool discard_when_basis_is_one
Definition: get_potential_grad.hpp:122
Definition: gradparx2.hpp:8
Definition: get_potential_grad.hpp:120
Kokkos::LayoutRight CLayout
Definition: space_settings.hpp:67
EMParField< DeviceType, DeviceOut > E_para_em
Definition: get_potential_grad.hpp:259
GradientMatrices< DeviceType > grad_matrices
Definition: get_potential_grad.hpp:71
GridField< DeviceOut, VT, PIT, TT, KT > gradient
Definition: get_potential_grad.hpp:253
bool potential_is_requested
Definition: get_potential_grad.hpp:256
Projection< HostType > half_plane_ff
Definition: grid.hpp:311
void get_field_grad(const Grid< DeviceType > &grid, const DomainDecomposition< DeviceType > &pol_decomp, const MagneticField< DeviceType > &magnetic_field, GetPotGradFieldArgs< DeviceIn, DeviceOut, VT, PIT, TT, KT > &args, GetPotentialGradTemp< DeviceType, DeviceOut > &tmp)
Definition: get_potential_grad.cpp:438
Definition: grid_field.hpp:25
View< double *, CLayout, Device, Kokkos::MemoryTraits< Kokkos::Unmanaged > > pot
Definition: get_potential_grad.hpp:125
AlternatingStorage< Device > gradient_z_alt
Definition: get_potential_grad.hpp:60
bool gradient_is_requested
Definition: get_potential_grad.hpp:257
View< double **, CLayout, DeviceIn, Kokkos::MemoryTraits< Kokkos::Unmanaged > > input_potential
Definition: get_potential_grad.hpp:252
View< double ***, CLayout, Device > gradient
Definition: get_potential_grad.hpp:62
PhiInterpType
Definition: globals.hpp:95
View< double ***, CLayout, Device > potential_rho
Definition: get_potential_grad.hpp:64
View< double *, CLayout, Device > right(int i_plane)
Definition: get_potential_grad.hpp:43
subroutine grid_deriv(grid, qty, qty_deriv_x, qty_deriv_y, psi_only)
Definition: search.F90:4063
constexpr PhiInterpType PIT_GLOBAL
Definition: globals.hpp:103
AlternatingStorage< Device > gradient_r_alt
Definition: get_potential_grad.hpp:59
double bt_sign
Whether toroidal field is reversed?
Definition: magnetic_field.hpp:19
View< double ***, CLayout, Device > field_rho
Definition: get_potential_grad.hpp:215
Field00< DeviceType > field00
Definition: get_potential_grad.hpp:258
GetPotGradFieldArgs(const View< double *, CLayout, DeviceIn > &input_potential_in, bool ignore_poloidal_dpot_in=false)
Definition: get_potential_grad.hpp:270
EMParField()
Definition: get_potential_grad.hpp:219
void request_potential(const GridField< DeviceOut, VarType::Scalar, PIT, TT, KT > &potential_in)
Definition: get_potential_grad.hpp:277
rtype * get_psn_pbd0_2_iseg_loc()
GyroAverageMatrices< DeviceType > gyro_avg_matrices
Definition: get_potential_grad.hpp:70
View< double *, CLayout, Device > z
Definition: get_potential_grad.hpp:127
Field00(const GridField< DeviceIn, VarType::Scalar, PhiInterpType::None, TorType::OnePlane, KinType::DriftKin > &field_in, bool discard_when_basis_is_one_in)
Definition: get_potential_grad.hpp:134
int nnode
Definition: get_potential_grad.hpp:52
Definition: field.hpp:50
void get_field_Ah_cv_ff(const Grid< DeviceType > &grid, const DomainDecomposition< DeviceType > &pol_decomp, GridField< HostType, VarType::Scalar, PhiInterpType::Planes, TorType::OnePlane, KinType::DriftKin > &Ah_cv_ff)
Definition: get_potential_grad.cpp:10
Definition: get_potential_grad.hpp:211
void request(const GridField< DeviceOut, VarType::Scalar, PhiInterpType::Planes, TorType::OnePlane, KinType::GyroKin > &output_field)
Definition: get_potential_grad.hpp:223
View< T *, CLayout, Device > my_mirror_view(const View< T *, CLayout, Device > &view, Device nd)
Definition: my_mirror_view.hpp:14
GridField< DeviceOut, VarType::Scalar, PhiInterpType::Planes, TorType::OnePlane, KinType::GyroKin > field_out
Definition: get_potential_grad.hpp:217
Kokkos::View< T *, Kokkos::LayoutRight, Device > my_subview(const Kokkos::View< T ****, Kokkos::LayoutRight, Device > &view, int i, int j, int k)
Definition: my_subview.hpp:8
Definition: magnetic_field.F90:1
Kokkos::Device< ExSpace, MemSpace > DeviceType
Definition: space_settings.hpp:47
GridField< DeviceOut, VarType::Scalar, PIT, TT, KT > potential
Definition: get_potential_grad.hpp:254
void calculate(const Grid< DeviceType > &grid, const View< double **, CLayout, Device > &field_para)
Definition: get_potential_grad.hpp:230
Definition: get_potential_grad.hpp:28
AlternatingStorage< Device > potential_alt
Definition: get_potential_grad.hpp:58
void request_gradient(const GridField< DeviceOut, VT, PIT, TT, KT > &gradient_in)
Definition: get_potential_grad.hpp:282
void parallel_for(const std::string name, int n_ptl, Function func, Option option, HostAoSoA aosoa_h, DeviceAoSoA aosoa_d)
Definition: streamed_parallel_for.hpp:252
void calculate_gradient(const Grid< DeviceType > &grid, const GradientMatrices< DeviceType > &grad_matrices)
Definition: get_potential_grad.hpp:146
View< double **, CLayout, Device > field
Definition: get_potential_grad.hpp:214
int nnode
Number of grid nodes.
Definition: grid.hpp:268
GradParXTmp gptx
Definition: get_potential_grad.hpp:73
View< double *, CLayout, Device > r
Definition: get_potential_grad.hpp:126
View< double **, CLayout, Device > field
Definition: get_potential_grad.hpp:29
void set_field(const Grid< DeviceType > &grid, const FieldFollowingCoordinates &ff, View< Field< VarType::Vector2D, PIT > *, CLayout, DeviceOut > &field00_ff_h)
Definition: get_potential_grad.hpp:162
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
bool requested
Definition: get_potential_grad.hpp:212
GetPotGradFieldArgs(const View< double **, CLayout, DeviceIn, Kokkos::MemoryTraits< Kokkos::Unmanaged >> &input_potential_in, bool ignore_poloidal_dpot_in=false)
Definition: get_potential_grad.hpp:262
View< double ****, CLayout, Device > gradient_rho
Definition: get_potential_grad.hpp:63
static int GPTLstop(const char *name)
Definition: timer_macro.hpp:10
Definition: get_potential_grad.hpp:50