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