7 extern "C" void my_qsort(
int left,
int right,
int*
idx,
double* psi_surf2_tmp,
int npsi_surf2);
9 extern "C" void generate_bicub_coeffs(
double* xv,
double* yv,
double* xc,
double* yc,
double* fval,
double* acoef );
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);
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);
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;
31 double det = 1.0/( dx1r*dx2z - dx2r*dx1z);
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;
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;
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 ){
84 int get_n_wall_nodes(
const View<int*, HostType>& rgn){
86 for(
int i = 0; i<rgn.size(); i++){
87 if(rgn[i]==Wall) count++;
92 printf(
"\ninit_wall: No wall nodes found!\n");
98 void setup_wall(
const View<int*, CLayout, HostType>& rgn, View<int*, CLayout, HostType>& wall_nodes, View<int*, CLayout, HostType>& node_to_wall){
100 for(
int i = 0; i<rgn.size(); i++){
102 wall_nodes[count]=i + 1;
103 node_to_wall[i]=count + 1;
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++){
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
124 && (rgn[i] != PrivateRegion || !exclude_private)
125 && (rgn[i] != Wall) && grad_psitheta ){
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);
139 Kokkos::View<double*,Kokkos::LayoutRight,Device> psi(
"psi", nnode);
142 magnetic_field.
psi_bicub.der_zero(gx(idx).r,gx(idx).z, psi_local);
143 psi(idx) = max(1e-70, psi_local);
157 template<
class Device>
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;
167 Kokkos::View<double*,Kokkos::LayoutRight,Device> psi_surf(
"psi_surf", nsurfaces);
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);
177 psi_surf(i) /= surf_vol;
183 Kokkos::View<int*,Kokkos::LayoutRight,Device> use_surf(
"use_surf", nsurfaces);
185 for(
int i=0; i<nsurfaces; i++){
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);
194 if (zmin <= use_surf_thresh){
200 psi_surf2_tmp = View<double*, HostType>(
NoInit(
"psi_surf2_tmp"),npsi_surf2);
202 for(
int i=0; i<nsurfaces; i++){
204 psi_surf2_tmp[j]=psi_surf(i);
210 View<int*, HostType>
idx(
NoInit(
"idx"),npsi_surf2);
211 for(
int i=0; i<npsi_surf2; i++) idx[i]=i;
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];
220 Kokkos::View<int**,Kokkos::LayoutRight,HostType> setup_guess_list_1d(
int npsi00,
double psi00min,
double dpsi00, View<double*, HostType>& psi_surf2){
222 Kokkos::View<int**,Kokkos::LayoutRight,HostType> guess_list_1d(
"guess_list_1d",2,npsi00);
227 for(
int i=0; i<npsi00; i++){
228 double psi = psi00min + i*dpsi00;
230 int up=psi_surf2.size();
231 int middle=down + (up-down)/2;
233 if (psi >= psi_surf2[middle-1]){
238 middle = down + (up-down)/2;
241 guess_list_1d(0,i)=down;
242 guess_list_1d(1,i)=up;
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;
250 for(
int i=0; i<(npsi00-1); i++){
254 guess_list_1d(1, i)=guess_list_1d(1, i+1);
257 return guess_list_1d;
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 ){
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);
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