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"
#include "natural_boundary.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;
    bool use_natural;
 
    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 && magnetic_field.is_below_xpt_tangent(r, z))
        || (grid.psi(i) < s.out_bd_psi_priv2 && magnetic_field.is_above_xpt2_tangent(r, z)))
       && !grid.node_is_in_region_1_or_2_no_wall(i) ) return true;
 
    // Use this condition only if the private region is also excluded
    if(s.excl_sep_leg && grid.midplane.node_is_on_separatrix_leg(magnetic_field, i)) 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);
                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:
 
    bool use_natural;
    NaturalBoundary natural; // Optional matrix setup for natural boundary settings
 
    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)),
          use_natural(settings.use_natural)
    {
        if(use_natural) natural = NaturalBoundary(grid, static_cast<VertexList>(*this));
    }
};
 
#endif