boundary.hpp

Go to the documentation of this file.

#ifndef BOUNDARY_HPP
#define BOUNDARY_HPP

#include "field.hpp"
#include "vertex_list.hpp"
#include "grid.hpp"

struct BoundarySettings{
    double in_bd_psi;
    double out_bd_psi;
    double out_bd_psi_priv1;
    double out_bd_psi_priv2;
    double bd_ext_near_wall;
    bool excl_inner;
    bool excl_outer;
    bool excl_private;
    bool excl_near_wall;
    bool excl_wall;
    bool excl_sep_leg;

    void extend(double inner_ext, double outer_ext){
        in_bd_psi += inner_ext;
        out_bd_psi -= outer_ext;
        out_bd_psi_priv1 += outer_ext;
        out_bd_psi_priv2 += outer_ext;
    }
};

static KOKKOS_INLINE_FUNCTION bool exclude_node(const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid,
                                         const View<int*,CLayout,DeviceType>& num_t_node, const View<int**,CLayout,DeviceType>& tr_node,
                                         const BoundarySettings& s, int i){
    constexpr double epsilon = 1.0e-8;
    if(s.excl_inner && grid.psi(i) - epsilon < s.in_bd_psi) return true;

    if(s.excl_outer && grid.psi(i) > s.out_bd_psi) return true;

    double r,z;
    grid.get_rz_coordinates(i, r, z);

    if((s.excl_private
        || (grid.psi(i) < s.out_bd_psi_priv1 && z < magnetic_field.equil.xpt.z )
        || (grid.psi(i) < s.out_bd_psi_priv2 && z > magnetic_field.equil.xpt2.z && magnetic_field.equil.set_xpt2) )
       && !grid.node_is_in_region_1_or_2_no_wall(i) ) return true;

    // Use this condition only if the private region is also excluded
    constexpr double epsilon2 = 1.0e-5;
    if(s.excl_sep_leg
       && (grid.node_is_in_region_1_or_2_no_wall(i) && grid.node_is_in_region_2_or_3_no_wall(i)) // is_in_region_2
       && (-(r-magnetic_field.equil.xpt.r)*magnetic_field.equil.xpt_slope + (z-magnetic_field.equil.xpt.z) < 0.0)
       && (grid.psi(i) < (magnetic_field.equil.xpt_psi*(1.0 + epsilon2)))) return true;

    if(s.excl_near_wall){
        // check if i-node is connected to wall
        for (int j=0; j<num_t_node(i); j++){
            int itr = tr_node(j,i) + 1;
            for(int ip=0; ip<3; ip++){
                int i_nd = grid.get_node_index(itr, ip);
                if(grid.midplane.node_is_on_wall(i_nd)) return true;
            }
        }
#ifndef XGCA
         // check if i-node is near wall (bd_ext_near_wall)
        if(s.bd_ext_near_wall>0.0){
            for(int j=0; j<grid.nwall; j++){
                RZPair x2 = grid.get_wall_rz(j);
                double dist = grid.midplane.get_dist2_from_node(i, x2.r, x2.z);
                if(dist < s.bd_ext_near_wall*s.bd_ext_near_wall){
                    return true;
                }
            }
        }
#endif
    } else if(s.excl_wall && grid.midplane.node_is_on_wall(i)) return true;

    return false;
}

inline VertexList get_excluded_vertex_list(const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, const View<int*,CLayout,DeviceType>& num_t_node, const View<int**,CLayout,DeviceType>& tr_node, const BoundarySettings& settings){
    // Create vertex list that excludes the boundary
    return VertexList(grid.nnode, KOKKOS_LAMBDA(const int i){
            return exclude_node(magnetic_field, grid, num_t_node, tr_node, settings, i);
        });
}

class Boundary : public VertexList{

    public:

    Boundary(){}

    Boundary(const MagneticField<DeviceType>& magnetic_field, const Grid<DeviceType>& grid, const View<int*,CLayout,DeviceType>& num_t_node, const View<int**,CLayout,DeviceType>& tr_node, const BoundarySettings& settings)
        : VertexList(get_excluded_vertex_list(magnetic_field, grid, num_t_node, tr_node, settings)){}
};

#endif