XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
varying_background.hpp
Go to the documentation of this file.
1 #ifndef VARYING_BACKGROUND_HPP
2 #define VARYING_BACKGROUND_HPP
3 
4 #include "magnetic_field.hpp"
5 #include "plasma.hpp"
6 
7 extern "C" double* get_col_vb_vol_loc();
8 
9 template<class Device>
11 
12  enum{
19  };
20 
21  enum{
25  N
26  };
27 
28  double pin;
29  double pout;
30  int m;
31  int mtheta;
32  double inv_dp;
33  double inv_dtheta;
34 
35  View<double**, CLayout, Device> vol;
36  View<double****, CLayout, Device> bg;
37 
38  public:
39 
40  int period;
41 
43 
45  nlr.use_namelist("col_param");
46  period = nlr.get<int>("col_vb_period", 1);
47  pin = nlr.get<double>("col_vb_pin", magnetic_field.inpsi);
48  pout = nlr.get<double>("col_vb_pout", std::min(magnetic_field.outpsi, 1.03));
49  m = nlr.get<int>("col_vb_m", 50);
50  mtheta = nlr.get<int>("col_vb_mtheta", 8);
51 
52  // Normalize
53  pin *= magnetic_field.equil.xpt_psi;
54  pout *= magnetic_field.equil.xpt_psi;
55 
56  double dp = (pout - pin)/m;
57  double dtheta = TWOPI/mtheta;
58  inv_dp = 1.0/dp;
59  inv_dtheta = 1.0/dtheta;
60 
61  // Allocate views
62  bg = View<double****, CLayout, Device>(NoInit("vb_bg"), n_nonadiabatic_species, m+1, mtheta, N);
63  vol = View<double**, CLayout, Device>(NoInit("vb_vol"), m+1, mtheta);
64  }
65 
66  KOKKOS_INLINE_FUNCTION bool is_in_range(const MagneticField<DeviceType>& magnetic_field, double r, double z, double psi) const {
67  return pin<psi && psi<pout && z>magnetic_field.equil.xpt_z && magnetic_field.equil.is_in_region_1_or_2(r,z,psi);
68  }
69 
70  KOKKOS_INLINE_FUNCTION void interp_bg_profile(double psi, double theta, int isp, double& den, double& temp, double& up) const {
71  double tmp = (psi-pin)*inv_dp;
72  int j = int(tmp);
73  j = min(m-1,max(0,j));
74  int l = ((int)(floor(theta*inv_dtheta + 0.5)))%mtheta;
75  double bb = tmp - j;
76  double aa = 1.0 - bb;
77 
78  den = bg(isp,j,l,Density)*aa + bg(isp,j+1,l,Density)*bb;
79  temp = bg(isp,j,l,Temp)*aa + bg(isp,j+1,l,Temp)*bb;
80 
81  j = floor((psi-pin)*inv_dp + 0.5);
82  up = bg(isp,j,l,Flow);
83  }
84 
85  // get ion/electron density and temperature for collision routine
87  View<double***, CLayout, DeviceType> dum(NoInit("dum"),m+1,mtheta,N);
88  Kokkos::deep_copy(bg, 1.0e-99); // to avoid divide by zero
89  Kokkos::deep_copy(dum, 1.0e-99); // to avoid divide by zero
90 
91  auto vb = *this; // Use a local shallow copy so that the lambda function knows to copy it
92 
93  species.for_all_particles("col_snapshot", KOKKOS_LAMBDA( const int idx ){
94  // Get index info
95  AoSoAIndices<Device> inds(idx);
96 
97  // Copy from AoSoA to SoA
98  SimdParticles part;
99  AoSoA_2_simd(part, species.ptl(), inds.a, inds.s);
100 
101  Simd<double> psi;
102  magnetic_field.get_psi(part.ph.x(),psi);
103 
104  Simd<double> Bmag;
105  magnetic_field.bmag_interpol(part.ph.v(), Bmag);
106 
107  Simd<double> theta;
108  magnetic_field.equil.get_theta(part.ph.x(), theta);
109 
110  FORCE_SIMD
111  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
112  if(part.gid[i_simd]>0){
113  if(vb.pin<psi[i_simd] && psi[i_simd]<vb.pout && magnetic_field.equil.is_in_region_1_or_2(part.ph.r[i_simd],part.ph.z[i_simd],psi[i_simd])){
114  //energy calculation
115  double weight = part.ct.w0[i_simd];
116  double upar = part.ph.rho[i_simd]*Bmag[i_simd]*species.c_m;
117  double energy = 0.5*species.mass*upar*upar + part.ct.mu[i_simd]*Bmag[i_simd];
118 
119  double tmp = (psi[i_simd] - vb.pin)*vb.inv_dp;
120  int j = (int)(tmp);
121  j = min(vb.m-1,max(0,j));
122  int l = ((int)(floor(theta[i_simd]*vb.inv_dtheta + 0.5)))%vb.mtheta;
123  double bb = tmp-j;
124  double aa = 1.0-bb;
125 
126  Kokkos::atomic_add(&(dum(j,l,Density1)), weight*aa);
127  Kokkos::atomic_add(&(dum(j,l,Energy)), weight*energy*aa);
128  Kokkos::atomic_add(&(dum(j,l,VPara)), weight*upar*aa);
129  Kokkos::atomic_add(&(dum(j+1,l,Density1)), weight*bb);
130  Kokkos::atomic_add(&(dum(j+1,l,Energy)), weight*energy*bb);
131  Kokkos::atomic_add(&(dum(j+1,l,VPara)), weight*upar*bb);
132 
133 
134  j = floor((psi[i_simd]-vb.pin)*vb.inv_dp + 0.5);
135  Kokkos::atomic_add(&(dum(j,l,VPara2)), weight);
136  Kokkos::atomic_add(&(dum(j,l,Density2)), weight*upar);
137  }
138  }
139  }
140  });
141 
142 #ifdef USE_MPI
143  TIMER("COL_SNAP_RED",
144  MPI_Allreduce(MPI_IN_PLACE, dum.data(), dum.size(), MPI_DOUBLE, MPI_SUM, SML_COMM_WORLD) );
145 #endif
146 
147  // Copy volume from fortran array to device
149 
150  int isp = species.nonadiabatic_idx;
151  int n = (m+1)*mtheta;
152  Kokkos::parallel_for("snapshot_convert", Kokkos::RangePolicy<ExSpace>(0, n), KOKKOS_LAMBDA( const int idx){
153  int i = idx%vb.mtheta; // Faster index
154  int j = idx/vb.mtheta;
155  constexpr double prefac = 2.0/3.0*J_2_EV; // Convert to eV
156  const double vpara_normed = dum(j,i,VPara)/dum(j,i,Density1);
157  const double epara = 0.5*species.mass*(vpara_normed*vpara_normed);
158  vb.bg(isp,j,i,Temp) = prefac*(dum(j,i,Energy) - epara)/dum(j,i,Density1);
159  vb.bg(isp,j,i,Density) = dum(j,i,Density1)/vb.vol(j,i);
160  vb.bg(isp,j,i,Flow) = dum(j,i,VPara2)/dum(j,i,Density2); // mean flow for moving frame collision
161  });
162  }
163 };
164 
165 #endif
void array_deep_copy(T *array, const Kokkos::View< T *, Kokkos::LayoutRight, Device > &view)
Definition: array_deep_copy.hpp:11
double inpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:25
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:353
Definition: varying_background.hpp:13
double c_m
c/m
Definition: species.hpp:77
VaryingBackground()
Definition: varying_background.hpp:42
View< double ****, CLayout, Device > bg
Density, Temperature and velocity of each shell for background update.
Definition: varying_background.hpp:36
Definition: varying_background.hpp:18
KOKKOS_INLINE_FUNCTION VecParticles * ptl() const
Definition: species.hpp:465
int m
Definition: varying_background.hpp:30
Definition: varying_background.hpp:14
Definition: varying_background.hpp:22
View< double **, CLayout, Device > vol
Volume.
Definition: varying_background.hpp:35
Definition: NamelistReader.hpp:163
Definition: magnetic_field.hpp:9
KOKKOS_INLINE_FUNCTION bool is_in_range(const MagneticField< DeviceType > &magnetic_field, double r, double z, double psi) const
Definition: varying_background.hpp:66
int a
The index in the inner array of the AoSoA.
Definition: particles.hpp:120
Definition: varying_background.hpp:15
KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector &v, Simd< double > &bmag) const
Definition: magnetic_field.tpp:185
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays) ...
Definition: species.hpp:72
int mtheta
Definition: varying_background.hpp:31
KOKKOS_INLINE_FUNCTION void AoSoA_2_simd(SimdParticles &part_one, const VecParticles *part, int a_vec, int s_vec)
Definition: particles.tpp:91
Simd< double > rho
Definition: particles.hpp:21
Definition: varying_background.hpp:24
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector2D &x, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:37
Definition: varying_background.hpp:16
double mass
Particle mass.
Definition: species.hpp:74
constexpr double J_2_EV
Conversion rate J to ev.
Definition: globals.hpp:120
void for_all_particles(const std::string label, F lambda_func) const
Definition: species.hpp:395
#define TIMER(N, F)
Definition: timer_macro.hpp:29
double outpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:26
idx
Definition: diag_f0_df_port1.hpp:32
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:322
Simd< double > r
Definition: particles.hpp:18
KOKKOS_INLINE_FUNCTION SimdVector2D & x()
Definition: particles.hpp:27
KOKKOS_INLINE_FUNCTION SimdVector & v()
Definition: particles.hpp:39
KOKKOS_INLINE_FUNCTION void interp_bg_profile(double psi, double theta, int isp, double &den, double &temp, double &up) const
Definition: varying_background.hpp:70
SimdPhase ph
Definition: particles.hpp:59
Definition: varying_background.hpp:17
double pin
Definition: varying_background.hpp:28
void update(const MagneticField< DeviceType > &magnetic_field, Species< DeviceType > &species)
Definition: varying_background.hpp:86
Definition: particles.hpp:58
int period
Period of background updating.
Definition: varying_background.hpp:40
int s
The index in the outer array of the AoSoA.
Definition: particles.hpp:119
#define FORCE_SIMD
Definition: simd.hpp:9
Simd< double > z
Definition: particles.hpp:19
constexpr double TWOPI
Definition: globals.hpp:123
Definition: varying_background.hpp:23
Simd< long long int > gid
Definition: particles.hpp:61
Equilibrium< Device > equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
Simd< double > w0
Definition: particles.hpp:53
Definition: magnetic_field.F90:1
SimdConstants ct
Definition: particles.hpp:60
Definition: species.hpp:66
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
double pout
Definition: varying_background.hpp:29
Simd< double > mu
Definition: particles.hpp:52
Definition: varying_background.hpp:10
double inv_dtheta
Definition: varying_background.hpp:33
double * get_col_vb_vol_loc()
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:61
Definition: particles.hpp:118
Definition: varying_background.hpp:25
double inv_dp
Definition: varying_background.hpp:32
VaryingBackground(NLReader::NamelistReader &nlr, const MagneticField< DeviceType > &magnetic_field, int n_nonadiabatic_species)
Definition: varying_background.hpp:44