XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
grid_setup.hpp
Go to the documentation of this file.
1 #ifndef GRID_SETUP_HPP
2 #define GRID_SETUP_HPP
3 
4 #include "globals.hpp"
5 #include "grid.hpp"
6 
7 extern "C" void my_qsort(int left,int right,int* idx,double* psi_surf2_tmp,int npsi_surf2);
8 
9 extern "C" void generate_bicub_coeffs(double* xv, double* yv, double* xc, double* yc, double* fval, double* acoef );
10 
11 
12 extern "C" void init_guess_xtable_fort( int* nguess_2, double* grid_guess_max, double* grid_guess_min, double* grid_inv_guess_d, int* ntriangle, int* nnodes, const RZPair* grid_x,
13  Vertex* grid_nd, int* guess_xtable, int* guess_count, int* guess_list_len_in);
14 
15 
16 extern "C" void init_guess_list_fort( int* nguess_2, double* grid_guess_max, double* grid_guess_min, double* grid_inv_guess_d, int* ntriangle, int* nnodes, const RZPair* grid_x,
17  Vertex* grid_nd, int* guess_list, int* guess_xtable, int* guess_count, int* guess_list_len_in);
18 
19 // Provides setup and checks for an analytic example of a grid
20 
21 namespace{
22 
23 // Mapping matrix init
24 void setup_mapping(const View<Vertex*, CLayout, HostType> nd, const View<RZPair*, CLayout, HostType> x, View<VertexMap*, CLayout, HostType>& mapping){
25  Kokkos::parallel_for("setup_mapping", Kokkos::RangePolicy<HostExSpace>(0,mapping.extent(0)), KOKKOS_LAMBDA( const int i ){
26  double dx1r=x[nd[i].vertex[0]-1].r - x[nd[i].vertex[2]-1].r;
27  double dx1z=x[nd[i].vertex[0]-1].z - x[nd[i].vertex[2]-1].z;
28  double dx2r=x[nd[i].vertex[1]-1].r - x[nd[i].vertex[2]-1].r;
29  double dx2z=x[nd[i].vertex[1]-1].z - x[nd[i].vertex[2]-1].z;
30 
31  double det = 1.0/( dx1r*dx2z - dx2r*dx1z);
32 
33  mapping[i].vertex[0].r = dx2z*det;
34  mapping[i].vertex[1].r =-dx2r*det;
35  mapping[i].vertex[0].z =-dx1z*det;
36  mapping[i].vertex[1].z = dx1r*det;
37 
38 
39  mapping[i].vertex[2].r = x[nd[i].vertex[2]-1].r;
40  mapping[i].vertex[2].z = x[nd[i].vertex[2]-1].z;
41  });
42 }
43 
44 // These are the regions designated in the .node file
46  Inside=0,
47  OnWall
48 };
49 
50 // These are the regions used in XGC
51 enum Region{
52  RegionOne=1,
53  RegionTwo,
54  PrivateRegion,
55  Wall=100 // temporary, since 100 is the value used elsewhere
56 };
57 
58 void set_psi(){
59  //for(int i = 0; i<psi.size(); i++){
60  //grid%psi(i)=psi_interpol(grid%x(1,i),grid%x(2,i),0,0);
61 }
62 
63 void set_rgn(View<int*, CLayout, HostType>& rgn){
64  Kokkos::parallel_for("set_rgn", Kokkos::RangePolicy<HostExSpace>(0,rgn.extent(0)), KOKKOS_LAMBDA( const int i ){
65  if(rgn[i]==Inside){
66  /*
67  if(is_rgn12(grid%x(1,i),grid%x(2,i),grid%psi(i)) .and. &
68  grid%psi(i) > eq_x_psi -epsil_psi ) then
69  grid%rgn(i)=2
70  elseif( is_rgn1(grid%x(1,i),grid%x(2,i),grid%psi(i)) ) then
71  grid%rgn(i)=1
72  else
73  grid%rgn(i)=3
74  endif
75  */
76  rgn[i]=RegionOne; // dummy
77  } else {
78  rgn[i]=Wall;
79  }
80  });
81 }
82 
83 // Counts number of wall nodes
84 int get_n_wall_nodes(const View<int*, HostType>& rgn){
85  int count=0;
86  for(int i = 0; i<rgn.size(); i++){
87  if(rgn[i]==Wall) count++;
88  }
89 
90  if (count == 0){
91  // No wall nodes found --> print warning and return
92  printf("\ninit_wall: No wall nodes found!\n");
93  }
94 
95  return count;
96 }
97 
98 void setup_wall(const View<int*, CLayout, HostType>& rgn, View<int*, CLayout, HostType>& wall_nodes, View<int*, CLayout, HostType>& node_to_wall){
99  int count=0;
100  for(int i = 0; i<rgn.size(); i++){
101  if(rgn[i]==Wall){
102  wall_nodes[count]=i + 1; // 1-indexed
103  node_to_wall[i]=count + 1; // 1-indexed
104  count++;
105  }else{
106  node_to_wall[i]=0;
107  }
108  }
109 }
110 
111 template<class Device>
112 View<int*,HostType> setup_basis(const MagneticField<Device>& magnetic_field, const View<RZPair*,HostType>& x, const View<int*,HostType>& rgn, double gradpsi_limit, bool exclude_private, bool grad_psitheta){
113  View<int*,HostType> basis(NoInit("basis"),rgn.size());
114  const double eps=1.0e-4;
115  for(int i = 0; i<basis.size(); i++){
116  double r = x[i].r;
117  double z = x[i].z;
118  double d1=sqrt((r - magnetic_field.equil.axis_r)*(r - magnetic_field.equil.axis_r) + (z - magnetic_field.equil.axis_z)*(z - magnetic_field.equil.axis_z));
119  double d2=sqrt((r - magnetic_field.equil.xpt_r) *(r - magnetic_field.equil.xpt_r) + (z - magnetic_field.equil.xpt_z) *(z - magnetic_field.equil.xpt_z) );
120  double dpsi_dr, dpsi_dz, psi_local;
121  magnetic_field.psi_bicub.der_one(r,z, psi_local,dpsi_dr,dpsi_dz );
122  double gradpsi = sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
123  if (d1 > eps && d2 > eps && gradpsi > gradpsi_limit // If not too close to the xpoint or the axis, and gradpsi is large enough
124  && (rgn[i] != PrivateRegion || !exclude_private) // If not in private region if that's excluded
125  && (rgn[i] != Wall) && grad_psitheta ){ // If we're not at the wall, and if grad_psitheta is true
126  basis[i]=0;
127  } else {
128  basis[i]=1;
129  }
130  }
131  return basis;
132 }
133 
134 template<class Device>
135 Kokkos::View<double*,Kokkos::LayoutRight,Device> setup_grid_psi(const MagneticField<Device>& magnetic_field, Kokkos::View<RZPair*,Kokkos::LayoutRight,Device>& gx){
136  int nnode = gx.extent(0);
137 
138  // Calculate psi at grid points
139  Kokkos::View<double*,Kokkos::LayoutRight,Device> psi("psi", nnode);
140  Kokkos::parallel_for("psi_calc", Kokkos::RangePolicy<ExSpace>(0, nnode), KOKKOS_LAMBDA( const int idx ){
141  double psi_local;
142  magnetic_field.psi_bicub.der_zero(gx(idx).r,gx(idx).z, psi_local);
143  psi(idx) = max(1e-70, psi_local);
144  });
145 
146  Kokkos::fence();
147 
148  return psi;
149 }
150 
151 /*
152  int nsurfaces;
153  std::vector<int> nnodes_on_surface;
154  int sum_nnodes_on_surfaces;
155  std::vector<int> nodes_of_surfaces;
156  */
157 template<class Device>
158 void setup_psi_surf2(const MagneticField<Device>& magnetic_field, int nsurfaces,
159  const View<int*, HostType>& nnodes_on_surface, const View<View<int*,HostType>*,HostType>& surf_idx,
160  const Kokkos::View<double*,Kokkos::LayoutRight,Device>& psi,
161  const Kokkos::View<double*,Kokkos::LayoutRight,HostType>& node_vol_h,
162  const Kokkos::View<RZPair*,Kokkos::LayoutRight,Device>& gx,
163  View<double*, HostType>& psi_surf2_tmp){
164  const double use_surf_thresh=3.0e-2;
165 
166  // Initializes to 0.0
167  Kokkos::View<double*,Kokkos::LayoutRight,Device> psi_surf("psi_surf", nsurfaces);
168 
169  // Calculate psi on the flux-surfaces
170  for(int i = 0; i<nsurfaces; i++){
171  double surf_vol = 0.0;
172  for(int j=0; j<nnodes_on_surface[i]; j++){
173  int k=surf_idx(i)(j);
174  psi_surf(i) += psi(k)*node_vol_h(k);
175  surf_vol += node_vol_h(k);
176  }
177  psi_surf(i) /= surf_vol;
178  }
179 
180 
181  // Determine which flux-surfaces to use for 1D diagnosis
182  // Using int instead of bool to initialize to 0
183  Kokkos::View<int*,Kokkos::LayoutRight,Device> use_surf("use_surf", nsurfaces);
184  int npsi_surf2=0;
185  for(int i=0; i<nsurfaces; i++){
186  // Test whether surface crosses outer midplane
187  double zmin=1e10;
188  for(int j=0; j<nnodes_on_surface[i]; j++){
189  int k=surf_idx(i)(j);
190  if ( (std::abs(gx(k).z-magnetic_field.equil.axis_z) < zmin) && (gx(k).r-(magnetic_field.equil.axis_r - 5.0e-5) >= 0.0) ){
191  zmin=std::abs(gx(k).z -magnetic_field.equil.axis_z);
192  }
193  }
194  if (zmin <= use_surf_thresh){
195  use_surf(i)=1;
196  npsi_surf2++;
197  }
198  }
199 
200  psi_surf2_tmp = View<double*, HostType>(NoInit("psi_surf2_tmp"),npsi_surf2);
201  int j=0;
202  for(int i=0; i<nsurfaces; i++){
203  if (use_surf(i)==1){
204  psi_surf2_tmp[j]=psi_surf(i);
205  j++;
206  }
207  }
208 
209  // Make sure surfaces are in ascending order
210  View<int*, HostType> idx(NoInit("idx"),npsi_surf2);
211  for(int i=0; i<npsi_surf2; i++) idx[i]=i;
212  int left=1;
213  int right=npsi_surf2;
214  my_qsort(left,right,idx.data(),psi_surf2_tmp.data(),npsi_surf2);
215  View<double*, HostType> tmp(NoInit("tmp"),npsi_surf2);
216  for(int i=0; i<npsi_surf2; i++) tmp[i]=psi_surf2_tmp[idx[i]];
217  for(int i=0; i<npsi_surf2; i++) psi_surf2_tmp[i]=tmp[i];
218 }
219 
220 Kokkos::View<int**,Kokkos::LayoutRight,HostType> setup_guess_list_1d(int npsi00, double psi00min, double dpsi00, View<double*, HostType>& psi_surf2){
221  // Allocate view
222  Kokkos::View<int**,Kokkos::LayoutRight,HostType> guess_list_1d("guess_list_1d",2,npsi00);
223 
224  // Run a binary search for each psi-value of
225  // the uniform psi00 grid --> upper and lower
226  // psi-index on the non-uniform mesh
227  for(int i=0; i<npsi00; i++){
228  double psi = psi00min + i*dpsi00;
229  int down=1;
230  int up=psi_surf2.size();
231  int middle=down + (up-down)/2;
232  while(up-down > 1){
233  if (psi >= psi_surf2[middle-1]){
234  down = middle;
235  }else{
236  up = middle;
237  }
238  middle = down + (up-down)/2;
239  }
240 
241  guess_list_1d(0,i)=down;
242  guess_list_1d(1,i)=up;
243 
244  if (psi < psi_surf2[down-1] || psi > psi_surf2[up-1]){
245  guess_list_1d(0,i) = -1;
246  guess_list_1d(1,i) = -1;
247  }
248  }
249 
250  for(int i=0; i<(npsi00-1); i++){
251  // Correct the upper bound to make sure that
252  // upper and lower bound enclose the interval
253  // psi00min+i*dpsi00 to psi00min+(i+1)*dpsi00
254  guess_list_1d(1, i)=guess_list_1d(1, i+1);
255  }
256 
257  return guess_list_1d;
258 }
259 
260 // Calculate and store magnetic field at nodes
261 void setup_grid_bfield(const MagneticField<DeviceType>& magnetic_field, const Kokkos::View<RZPair*,Kokkos::LayoutRight,DeviceType>& gx, Kokkos::View<double*,Kokkos::LayoutRight,DeviceType>& bfield){
262  Kokkos::parallel_for("node_bfield", Kokkos::RangePolicy<ExSpace>(0,bfield.extent(0)), KOKKOS_LAMBDA( const int idx ){
263  double br, bz, bphi;
264  magnetic_field.bvec_interpol(gx(idx).r,gx(idx).z,0.0,br,bz,bphi);
265  bfield(idx) = sqrt(br*br + bz*bz + bphi*bphi);
266  });
267 }
268 
269 /*void plot_grid(Grid<HostType>& grid){
270  int nx = 80;
271  int ny = 40;
272  double dr = (equil.max_r - equil.min_r)/nx;
273  double dz = (equil.max_z - equil.min_z)/ny;
274  printf("\n");
275  std::vector<int> prevrow(nx+1,-1);
276 
277  for (int j=ny;j>=0;j--){
278  printf("\n");
279  int prevnum = -1;
280  for (int i=0;i<=nx;i++){
281  double r = equil.min_r + i*dr;
282  double z = equil.min_z + j*dz;
283  int num = get_analytic_psi(equil,r,z)*5.0;
284  char mychar = 48+num;
285  if (prevnum!=num || prevrow[i]!=num) printf("%c",mychar);
286  else printf(" ");
287  prevnum = num;
288  prevrow[i] = num;
289  }
290  }
291  printf("\n");
292 }
293 */
294 
295 }
296 #endif
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:150
Kokkos::Device< HostExSpace, HostMemSpace > HostType
Definition: space_settings.hpp:56
void init_guess_xtable_fort(int *nguess_2, double *grid_guess_max, double *grid_guess_min, double *grid_inv_guess_d, int *ntriangle, int *nnodes, const RZPair *grid_x, Vertex *grid_nd, int *guess_xtable, int *guess_count, int *guess_list_len_in)
void my_qsort(int left, int right, int *idx, double *psi_surf2_tmp, int npsi_surf2)
Definition: magnetic_field.hpp:9
Definition: grid_structs.hpp:28
idx
Definition: diag_f0_df_port1.hpp:32
FileRegion
Definition: grid_setup.hpp:45
void generate_bicub_coeffs(double *xv, double *yv, double *xc, double *yc, double *fval, double *acoef)
void init_guess_list_fort(int *nguess_2, double *grid_guess_max, double *grid_guess_min, double *grid_inv_guess_d, int *ntriangle, int *nnodes, const RZPair *grid_x, Vertex *grid_nd, int *guess_list, int *guess_xtable, int *guess_count, int *guess_list_len_in)
Definition: grid_structs.hpp:24
Region
Definition: grid_setup.hpp:51
Equilibrium< Device > equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
Definition: magnetic_field.F90:1
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
Bicub< Device > psi_bicub
The object for interpolating psi (magnetic flux surfaces)
Definition: magnetic_field.hpp:30
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:61