grid_setup.hpp

Go to the documentation of this file.

#ifndef GRID_SETUP_HPP
#define GRID_SETUP_HPP

#include "globals.hpp"
#include "grid.hpp"

extern "C" void my_qsort(int left,int right,int* idx,double* psi_surf2_tmp,int npsi_surf2);

extern "C" void generate_bicub_coeffs(double* xv, double* yv, double* xc, double* yc, double* fval, double* acoef );


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,
                            Vertex* grid_nd, int* guess_xtable, int* guess_count, int* guess_list_len_in);


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,
                                     Vertex* grid_nd, int* guess_list, int* guess_xtable, int* guess_count, int* guess_list_len_in);

// Provides setup and checks for an analytic example of a grid

namespace{

// Mapping matrix init
void setup_mapping(const View<Vertex*, CLayout, HostType> nd, const View<RZPair*, CLayout, HostType> x, View<VertexMap*, CLayout, HostType>& mapping){
    Kokkos::parallel_for("setup_mapping", Kokkos::RangePolicy<HostExSpace>(0,mapping.extent(0)), KOKKOS_LAMBDA( const int i ){
       double dx1r=x[nd[i].vertex[0]-1].r - x[nd[i].vertex[2]-1].r;
       double dx1z=x[nd[i].vertex[0]-1].z - x[nd[i].vertex[2]-1].z;
       double dx2r=x[nd[i].vertex[1]-1].r - x[nd[i].vertex[2]-1].r;
       double dx2z=x[nd[i].vertex[1]-1].z - x[nd[i].vertex[2]-1].z;

       double det = 1.0/( dx1r*dx2z - dx2r*dx1z);

       mapping[i].vertex[0].r = dx2z*det;
       mapping[i].vertex[1].r =-dx2r*det;
       mapping[i].vertex[0].z =-dx1z*det;
       mapping[i].vertex[1].z = dx1r*det;


       mapping[i].vertex[2].r = x[nd[i].vertex[2]-1].r;
       mapping[i].vertex[2].z = x[nd[i].vertex[2]-1].z;
    });
}

// These are the regions designated in the .node file
enum FileRegion{
    Inside=0,
    OnWall
};

// These are the regions used in XGC
enum Region{
    RegionOne=1,
    RegionTwo,
    PrivateRegion,
    Wall=100 // temporary, since 100 is the value used elsewhere
};

void set_psi(){
    //for(int i = 0; i<psi.size(); i++){
    //grid%psi(i)=psi_interpol(grid%x(1,i),grid%x(2,i),0,0);
}

void set_rgn(View<int*, CLayout, HostType>& rgn){
    Kokkos::parallel_for("set_rgn", Kokkos::RangePolicy<HostExSpace>(0,rgn.extent(0)), KOKKOS_LAMBDA( const int i ){
        if(rgn[i]==Inside){
            /*
            if(is_rgn12(grid%x(1,i),grid%x(2,i),grid%psi(i)) .and. &
                grid%psi(i) > eq_x_psi -epsil_psi ) then
                grid%rgn(i)=2
             elseif( is_rgn1(grid%x(1,i),grid%x(2,i),grid%psi(i)) ) then
                grid%rgn(i)=1
             else
                grid%rgn(i)=3
             endif
             */
             rgn[i]=RegionOne; // dummy
        } else {
            rgn[i]=Wall;
        }
    });
}

// Counts number of wall nodes
int get_n_wall_nodes(const View<int*, HostType>& rgn){
    int count=0;
    for(int i = 0; i<rgn.size(); i++){
        if(rgn[i]==Wall) count++;
    }

    if (count == 0){
      // No wall nodes found --> print warning and return
      printf("\ninit_wall: No wall nodes found!\n");
    }

    return count;
}

void setup_wall(const View<int*, CLayout, HostType>& rgn, View<int*, CLayout, HostType>& wall_nodes, View<int*, CLayout, HostType>& node_to_wall){
    int count=0;
    for(int i = 0; i<rgn.size(); i++){
        if(rgn[i]==Wall){
            wall_nodes[count]=i + 1; // 1-indexed
            node_to_wall[i]=count + 1; // 1-indexed
            count++;
        }else{
            node_to_wall[i]=0;
        }
    }
}

template<class Device>
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){
    View<int*,HostType> basis(NoInit("basis"),rgn.size());
    const double eps=1.0e-4;
    for(int i = 0; i<basis.size(); i++){
        double r = x[i].r;
        double z = x[i].z;
        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));
        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) );
        double dpsi_dr, dpsi_dz, psi_local;
        magnetic_field.psi_bicub.der_one(r,z, psi_local,dpsi_dr,dpsi_dz );
        double gradpsi = sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
        if (d1 > eps && d2 > eps && gradpsi > gradpsi_limit // If not too close to the xpoint or the axis, and gradpsi is large enough
          && (rgn[i] != PrivateRegion || !exclude_private) // If not in private region if that's excluded
          && (rgn[i] != Wall) && grad_psitheta ){ // If we're not at the wall, and if grad_psitheta is true
            basis[i]=0;
        } else {
            basis[i]=1;
        }
    }
    return basis;
}

template<class Device>
Kokkos::View<double*,Kokkos::LayoutRight,Device> setup_grid_psi(const MagneticField<Device>& magnetic_field, Kokkos::View<RZPair*,Kokkos::LayoutRight,Device>& gx){
    int nnode = gx.extent(0);

    // Calculate psi at grid points
    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi("psi", nnode);
    Kokkos::parallel_for("psi_calc", Kokkos::RangePolicy<ExSpace>(0, nnode), KOKKOS_LAMBDA( const int idx ){
        double psi_local;
        magnetic_field.psi_bicub.der_zero(gx(idx).r,gx(idx).z, psi_local);
        psi(idx) = max(1e-70, psi_local);
    });

    Kokkos::fence();

    return psi;
}

/*
 int nsurfaces;
    std::vector<int> nnodes_on_surface;
    int sum_nnodes_on_surfaces;
    std::vector<int> nodes_of_surfaces;
    */
template<class Device>
void setup_psi_surf2(const MagneticField<Device>& magnetic_field, int nsurfaces,
                     const View<int*, HostType>& nnodes_on_surface, const View<View<int*,HostType>*,HostType>& surf_idx,
                     const Kokkos::View<double*,Kokkos::LayoutRight,Device>& psi,
                     const Kokkos::View<double*,Kokkos::LayoutRight,HostType>& node_vol_h,
                     const Kokkos::View<RZPair*,Kokkos::LayoutRight,Device>& gx, 
                     View<double*, HostType>& psi_surf2_tmp){
    const double use_surf_thresh=3.0e-2;

    // Initializes to 0.0
    Kokkos::View<double*,Kokkos::LayoutRight,Device> psi_surf("psi_surf", nsurfaces);

    // Calculate psi on the flux-surfaces
    for(int i = 0; i<nsurfaces; i++){
        double surf_vol = 0.0;
        for(int j=0; j<nnodes_on_surface[i]; j++){
            int k=surf_idx(i)(j);
            psi_surf(i) += psi(k)*node_vol_h(k);
            surf_vol += node_vol_h(k);
        }
        psi_surf(i) /= surf_vol;
    }


    // Determine which flux-surfaces to use for 1D diagnosis
    // Using int instead of bool to initialize to 0
    Kokkos::View<int*,Kokkos::LayoutRight,Device> use_surf("use_surf", nsurfaces);
    int npsi_surf2=0;
    for(int i=0; i<nsurfaces; i++){
        // Test whether surface crosses outer midplane
        double zmin=1e10;
        for(int j=0; j<nnodes_on_surface[i]; j++){
            int k=surf_idx(i)(j);
            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) ){
                zmin=std::abs(gx(k).z -magnetic_field.equil.axis_z);
            }
        }
        if (zmin <= use_surf_thresh){
            use_surf(i)=1;
            npsi_surf2++;
        }
    }

    psi_surf2_tmp = View<double*, HostType>(NoInit("psi_surf2_tmp"),npsi_surf2);
    int j=0;
    for(int i=0; i<nsurfaces; i++){
        if (use_surf(i)==1){
            psi_surf2_tmp[j]=psi_surf(i);
            j++;
        }
    }

    // Make sure surfaces are in ascending order
    View<int*, HostType> idx(NoInit("idx"),npsi_surf2);
    for(int i=0; i<npsi_surf2; i++) idx[i]=i;
    int left=1;
    int right=npsi_surf2;
    my_qsort(left,right,idx.data(),psi_surf2_tmp.data(),npsi_surf2);
    View<double*, HostType> tmp(NoInit("tmp"),npsi_surf2);
    for(int i=0; i<npsi_surf2; i++) tmp[i]=psi_surf2_tmp[idx[i]];
    for(int i=0; i<npsi_surf2; i++) psi_surf2_tmp[i]=tmp[i];
}

Kokkos::View<int**,Kokkos::LayoutRight,HostType> setup_guess_list_1d(int npsi00, double psi00min, double dpsi00, View<double*, HostType>& psi_surf2){
    // Allocate view
    Kokkos::View<int**,Kokkos::LayoutRight,HostType> guess_list_1d("guess_list_1d",2,npsi00);

    // Run a binary search for each psi-value of
    // the uniform psi00 grid --> upper and lower
    // psi-index on the non-uniform mesh
    for(int i=0; i<npsi00; i++){
        double psi = psi00min + i*dpsi00;
        int down=1;
        int up=psi_surf2.size();
        int middle=down + (up-down)/2;
        while(up-down > 1){
            if (psi >= psi_surf2[middle-1]){
                down = middle;
            }else{
                up = middle;
            }
            middle = down + (up-down)/2;
        }

        guess_list_1d(0,i)=down;
        guess_list_1d(1,i)=up;

        if (psi < psi_surf2[down-1] || psi > psi_surf2[up-1]){
            guess_list_1d(0,i) = -1;
            guess_list_1d(1,i) = -1;
        }
    }

    for(int i=0; i<(npsi00-1); i++){
        // Correct the upper bound to make sure that
        // upper and lower bound enclose the interval
        // psi00min+i*dpsi00 to psi00min+(i+1)*dpsi00
        guess_list_1d(1, i)=guess_list_1d(1, i+1);
    }

    return guess_list_1d;
}

// Calculate and store magnetic field at nodes
void setup_grid_bfield(const MagneticField<DeviceType>& magnetic_field, const Kokkos::View<RZPair*,Kokkos::LayoutRight,DeviceType>& gx, Kokkos::View<double*,Kokkos::LayoutRight,DeviceType>& bfield){
    Kokkos::parallel_for("node_bfield", Kokkos::RangePolicy<ExSpace>(0,bfield.extent(0)), KOKKOS_LAMBDA( const int idx ){
        double br, bz, bphi;
        magnetic_field.bvec_interpol(gx(idx).r,gx(idx).z,0.0,br,bz,bphi);
        bfield(idx) = sqrt(br*br + bz*bz + bphi*bphi);
    });
}

/*void plot_grid(Grid<HostType>& grid){
    int nx = 80;
    int ny = 40;
    double dr = (equil.max_r - equil.min_r)/nx;
    double dz = (equil.max_z - equil.min_z)/ny;
    printf("\n");
    std::vector<int> prevrow(nx+1,-1);

    for (int j=ny;j>=0;j--){
        printf("\n");
        int prevnum = -1;
        for (int i=0;i<=nx;i++){
            double r = equil.min_r + i*dr;
            double z = equil.min_z + j*dz;
            int num = get_analytic_psi(equil,r,z)*5.0;
            char mychar = 48+num;
            if (prevnum!=num || prevrow[i]!=num) printf("%c",mychar);
            else printf(" ");
            prevnum = num;
            prevrow[i] = num;
        }
    }
    printf("\n");
}
*/

}
#endif