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