7 extern "C" void my_qsort(
int left,
int right,
int*
idx,
double* psi_surf2_tmp,
int npsi_surf2);
9 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,
10 Vertex* grid_nd,
int* guess_xtable,
int* guess_count,
int* guess_list_len_in);
13 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,
14 Vertex* grid_nd,
int* guess_list,
int* guess_xtable,
int* guess_count,
int* guess_list_len_in);
21 void setup_mapping(
const View<Vertex*, CLayout, HostType> nd,
const View<RZPair*, CLayout, HostType> x, View<VertexMap*, CLayout, HostType>& mapping){
22 Kokkos::parallel_for(
"setup_mapping", Kokkos::RangePolicy<HostExSpace>(0,mapping.extent(0)), KOKKOS_LAMBDA(
const int i ){
23 double dx1r=x[nd[i].vertex[0]-1].r - x[nd[i].vertex[2]-1].r;
24 double dx1z=x[nd[i].vertex[0]-1].z - x[nd[i].vertex[2]-1].z;
25 double dx2r=x[nd[i].vertex[1]-1].r - x[nd[i].vertex[2]-1].r;
26 double dx2z=x[nd[i].vertex[1]-1].z - x[nd[i].vertex[2]-1].z;
28 double det = 1.0/( dx1r*dx2z - dx2r*dx1z);
30 mapping[i].vertex[0].r = dx2z*det;
31 mapping[i].vertex[1].r =-dx2r*det;
32 mapping[i].vertex[0].z =-dx1z*det;
33 mapping[i].vertex[1].z = dx1r*det;
36 mapping[i].vertex[2].r = x[nd[i].vertex[2]-1].r;
37 mapping[i].vertex[2].z = x[nd[i].vertex[2]-1].z;
60 template<
class Device,
class Gr
idDevice>
62 const View<double*,CLayout,GridDevice>& psi_h, View<int*, CLayout, GridDevice>& rgn_h){
64 int nnode = gx_h.extent(0);
65 View<RZPair*,CLayout,Device> gx(
"gx",nnode);
66 View<double*,CLayout,Device> psi(
"psi",nnode);
67 View<int*,CLayout,Device> rgn(
"rgn",nnode);
68 Kokkos::deep_copy(gx, gx_h);
69 Kokkos::deep_copy(psi, psi_h);
70 Kokkos::deep_copy(rgn, rgn_h);
72 Kokkos::parallel_for(
"set_rgn", Kokkos::RangePolicy<typename Device::execution_space>(0,rgn.extent(0)), KOKKOS_LAMBDA(
const int i ){
89 Kokkos::deep_copy(gx_h, gx);
90 Kokkos::deep_copy(psi_h, psi);
91 Kokkos::deep_copy(rgn_h, rgn);
95 int get_n_wall_nodes(
const View<int*, HostType>& rgn){
97 for(
int i = 0; i<rgn.size(); i++){
98 if(rgn[i]==Wall) count++;
103 printf(
"\ninit_wall: No wall nodes found!\n");
109 void setup_wall(
const View<int*, CLayout, HostType>& rgn, View<int*, CLayout, HostType>& wall_nodes, View<int*, CLayout, HostType>& node_to_wall){
111 for(
int i = 0; i<rgn.size(); i++){
113 wall_nodes[count]=i + 1;
114 node_to_wall[i]=count + 1;
123 template<
class Device,
class Gr
idDevice>
125 int nnode = gx_h.extent(0);
126 View<RZPair*,CLayout,Device> gx(
"gx",nnode);
127 Kokkos::deep_copy(gx, gx_h);
130 View<double*,CLayout,Device> psi(
"psi", nnode);
131 Kokkos::parallel_for(
"psi_calc", Kokkos::RangePolicy<typename Device::execution_space>(0, nnode), KOKKOS_LAMBDA(
const int idx ){
133 magnetic_field.
psi_bicub.der_zero(gx(idx).r,gx(idx).z, psi_local);
134 psi(idx) = max(1e-70, psi_local);
140 double psi_zero_eps = 5.0e-3;
141 Kokkos::parallel_for(
"psi_calc_check_first", Kokkos::RangePolicy<typename Device::execution_space>(0, 1), KOKKOS_LAMBDA(
const int idx ){
142 double dr = gx(idx).r - magnetic_field.
equil.
axis_r;
143 double dz = gx(idx).z - magnetic_field.
equil.
axis_z;
144 if(sqrt(dr*dr + dz*dz) < psi_zero_eps){
151 View<double*,CLayout,GridDevice> psi_h(
"psi", nnode);
152 Kokkos::deep_copy(psi_h, psi);
158 template<
class Device>
160 const View<double*,CLayout,HostType>& psi_h,
161 const View<int*, CLayout, HostType>& nnodes_on_surface_h,
162 const View<int**,CLayout,HostType>& surf_idx_h,
163 const View<double*,CLayout,HostType>& node_vol_h){
166 View<double*,CLayout,Device> psi(
NoInit(
"psi"), psi_h.layout());
167 View<int*, CLayout,Device> nnodes_on_surface(
NoInit(
"nnodes_on_surface"),nnodes_on_surface_h.layout());
168 View<int**,CLayout,Device> surf_idx(
NoInit(
"surf_idx"),surf_idx_h.layout());
169 View<double*,CLayout,Device> node_vol(
NoInit(
"node_vol"),node_vol_h.layout());
170 Kokkos::deep_copy(psi, psi_h);
171 Kokkos::deep_copy(nnodes_on_surface, nnodes_on_surface_h);
172 Kokkos::deep_copy(node_vol, node_vol_h);
173 Kokkos::deep_copy(surf_idx, surf_idx_h);
176 View<double*,CLayout,Device> psi_surf(
NoInit(
"psi_surf"), nnodes_on_surface.size());
179 Kokkos::parallel_for(
"get_psi_surf", Kokkos::RangePolicy<typename Device::execution_space>(0,psi_surf.size()), KOKKOS_LAMBDA(
const int i ){
180 double surf_vol = 0.0;
181 double psi_sum = 0.0;
182 for(
int j=0; j<nnodes_on_surface(i); j++){
183 int k=surf_idx(i, j)-1;
187 psi_sum += psi(k)*node_vol(k);
188 surf_vol += node_vol(k);
190 psi_surf(i) = psi_sum/surf_vol;
194 View<double*,CLayout,HostType> psi_surf_h(
NoInit(
"psi_surf"), psi_surf.layout());
195 Kokkos::deep_copy(psi_surf_h, psi_surf);
205 template<
class Device>
207 const View<int*, CLayout, HostType>& nnodes_on_surface_h,
208 const View<int**,CLayout, HostType>& surf_idx_h,
209 const View<double*,CLayout,HostType>& psi_surf_h,
210 const View<RZPair*,CLayout,HostType>& gx_h){
211 const double use_surf_thresh=3.0e-2;
212 int nnode = gx_h.extent(0);
213 View<RZPair*,CLayout,Device> gx(
"gx",nnode);
214 View<int*, CLayout, Device> nnodes_on_surface(
NoInit(
"nnodes_on_surface"),nnodes_on_surface_h.layout());
215 View<int**,CLayout, Device> surf_idx(
NoInit(
"surf_idx"),surf_idx_h.layout());
216 Kokkos::deep_copy(nnodes_on_surface, nnodes_on_surface_h);
217 Kokkos::deep_copy(surf_idx, surf_idx_h);
218 Kokkos::deep_copy(gx, gx_h);
220 int nsurfaces = nnodes_on_surface.size();
223 View<bool*,CLayout,Device> use_surf(
"use_surf", nsurfaces);
224 Kokkos::parallel_reduce(
"get_psi_surf2", Kokkos::RangePolicy<typename Device::execution_space>(0,psi_surf_h.size()), KOKKOS_LAMBDA(
const int i,
int& npsi_surf2_l){
226 for(
int j=0; j<nnodes_on_surface(i); j++){
227 int k=surf_idx(i, j)-1;
230 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) ){
231 zmin=std::abs(gx(k).z -magnetic_field.
equil.
axis_z);
235 if (zmin <= use_surf_thresh){
245 View<bool*,CLayout,HostType> use_surf_h(
"use_surf", nsurfaces);
246 Kokkos::deep_copy(use_surf_h, use_surf);
249 View<double*, CLayout, HostType> psi_surf2_h(
NoInit(
"psi_surf2"),npsi_surf2);
251 for(
int i=0; i<nsurfaces; i++){
253 psi_surf2_h(j) = psi_surf_h(i);
259 View<int*, HostType>
idx(
NoInit(
"idx"),npsi_surf2);
260 for(
int i=0; i<npsi_surf2; i++)
idx(i) = i + 1;
262 int right=npsi_surf2;
263 my_qsort(left,right,idx.data(),psi_surf2_h.data(),npsi_surf2);
264 View<double*, HostType> tmp(
NoInit(
"tmp"),npsi_surf2);
265 for(
int i=0; i<npsi_surf2; i++) tmp(i)=psi_surf2_h(
idx(i) - 1);
266 for(
int i=0; i<npsi_surf2; i++) psi_surf2_h(i)=tmp(i);
272 template<
class Device,
class Gr
idDevice>
273 View<double*,CLayout,GridDevice> setup_grid_bfield(
const MagneticField<Device>& magnetic_field,
const View<RZPair*,CLayout,GridDevice>& gx_h){
274 int nnode = gx_h.extent(0);
275 View<RZPair*,CLayout,Device> gx(
"gx",nnode);
276 Kokkos::deep_copy(gx, gx_h);
279 View<double*,CLayout,Device> bfield(
"bfield", nnode);
280 Kokkos::parallel_for(
"node_bfield", Kokkos::RangePolicy<typename Device::execution_space>(0,nnode), KOKKOS_LAMBDA(
const int idx ){
282 magnetic_field.
bvec_interpol(gx(idx).r,gx(idx).z,0.0,br,bz,bphi);
283 bfield(idx) = sqrt(br*br + bz*bz + bphi*bphi);
287 View<double*,CLayout,GridDevice> bfield_h(
"bfield", nnode);
288 Kokkos::deep_copy(bfield_h, bfield);
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:222
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:12
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
double axis_z
z coordinate of axis
Definition: equil.hpp:95
Definition: grid_structs.hpp:28
idx
Definition: diag_f0_df_port1.hpp:32
FileRegion
Definition: grid_setup.hpp:42
double axis_r
r coordinate of axis
Definition: equil.hpp:94
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:48
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
double epsil_psi
Not sure?
Definition: equil.hpp:85
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r, double z, double psi) const
Definition: equil.tpp:102
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:84