Upwind_gradient.hpp

Go to the documentation of this file.

//
// Created by jstark on 17.05.21.
//
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_UPWIND_GRADIENT_HPP
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_UPWIND_GRADIENT_HPP

// Include standard libraries
#include <cmath>

// Include OpenFPM header files
#include "Grid/grid_dist_id.hpp"
#include "FD_simple.hpp"
#include "Eno_Weno.hpp"

template <typename field_type>
static field_type upwinding(field_type dplus, field_type dminus, int sign)
{
    field_type grad_upwind = 0.0;
    if (dplus * sign < 0
            && (dminus + dplus) * sign < 0) grad_upwind = dplus;
    
    else if (dminus * sign > 0
            && (dminus + dplus) * sign > 0) grad_upwind = dminus;
    
    else if (dminus * sign < 0
            && dplus * sign > 0) grad_upwind = 0;
    return grad_upwind;
}


template <size_t Field, size_t Sign, typename gridtype, typename keytype>
auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
{
//  std::cout << boost::typeindex::type_id_with_cvr<std::decay_t<decltype(first_node)>>() << std::endl;
    
    typedef typename std::decay_t<decltype(grid.template get<Field>(key))> field_type;
    field_type dplus = 0.0, dminus = 0.0;
    int sign_phi0 = grid.template get<Sign> (key);
    
    switch(order)
    {
        case 1:
            dplus  = FD_forward<Field>(grid, key, d);
            dminus = FD_backward<Field>(grid, key, d);
            break;
        case 3:
            dplus  = ENO_3_Plus<Field>(grid, key, d);
            dminus = ENO_3_Minus<Field>(grid, key, d);
            break;
        case 5:
            dplus  = WENO_5_Plus<Field>(grid, key, d);
            dminus = WENO_5_Minus<Field>(grid, key, d);
            break;
        default:
            auto &v_cl = create_vcluster();
            if (v_cl.rank() == 0) std::cout << "Order of accuracy chosen not valid. Using default order 1." <<
                        std::endl;
            dplus  = FD_forward<Field>(grid, key, d);
            dminus = FD_backward<Field>(grid, key, d);
            break;
    }
    
    return upwinding(dplus, dminus, sign_phi0);
}


// Use upwinding for inner grid points and one sided backward / forward stencil at border (if one_sided_BC=true)
template <size_t Field, size_t Sign, size_t Gradient, typename gridtype>
void upwind_gradient(gridtype & grid, const bool one_sided_BC, size_t order)
{
    grid.template ghost_get<Field>(KEEP_PROPERTIES);
    grid.template ghost_get<Sign>(KEEP_PROPERTIES);
    
    auto dom = grid.getDomainIterator();
    
    if (one_sided_BC)
    {
        while (dom.isNext())
        {
            auto key = dom.get();
            auto key_g = grid.getGKey(key);
            
            for(size_t d = 0; d < gridtype::dims; d++ )
            {
                // Grid nodes inside and distance from boundary > stencil-width
                if (key_g.get(d) > 2 && key_g.get(d) < grid.size(d) - 3) // if point lays with min. 3 nodes distance to
                    // boundary
                {
                    grid.template get<Gradient> (key) [d] = FD_upwind<Field, Sign>(grid, key, d, order);
                }
                    
                    // Grid nodes in stencil-wide vicinity to boundary
                else if (key_g.get(d) > 0 && key_g.get(d) < grid.size(d) - 1) // if point lays not on the grid boundary
                {
                    grid.template get<Gradient> (key) [d] = FD_upwind<Field, Sign>(grid, key, d, 1);
                }
                
                else if (key_g.get(d) == 0) // if point lays at left boundary, use right sided kernel
                {
                    grid.template get<Gradient> (key) [d] = FD_forward<Field>(grid, key, d);
                }
                
                else if (key_g.get(d) >= grid.size(d) - 1) // if point lays at right boundary, use left sided kernel
                {
                    grid.template get<Gradient> (key) [d] = FD_backward<Field>(grid, key, d);
                }
            }
            ++dom;
        }
    }
    
    else
    {
        while (dom.isNext())
        {
            auto key = dom.get();
            for(size_t d = 0; d < gridtype::dims; d++)
            {
                grid.template get<Gradient> (key) [d] = FD_upwind<Field, Sign>(grid, key, d, order);
            }
            ++dom;
        }
    }
}

template <typename gridtype>
static bool ghost_width_is_sufficient(gridtype & grid, size_t required_width)
{
    auto ghost = grid.getDecomposition().getGhost();
    int np_ghost[gridtype::dims];
    
    for (int d = 0 ; d < gridtype::dims ; d++)
    {
        np_ghost[d] = ghost.getHigh(d) / grid.getSpacing()[d];
    }
    
    int min_width = np_ghost[0];
    for (int d = 0 ; d < gridtype::dims ; d++)
    {
        if(np_ghost[d] < min_width) min_width = np_ghost[d];
    }
    
    return (min_width >= required_width);
}

template <size_t Field_in, size_t Sign, size_t Gradient_out, typename gridtype>
void get_upwind_gradient(gridtype & grid, const size_t order=5, const bool one_sided_BC=true)
{
    grid.template ghost_get<Field_in>();
    grid.template ghost_get<Sign>();
    
    if (!one_sided_BC)
    {
        auto &v_cl = create_vcluster();
        if (v_cl.rank() == 0)
        {
            // Check if ghost layer thick enough for stencil
            size_t stencil_width;
            (order > 1) ? stencil_width = 3 : stencil_width = 1;
            if (!ghost_width_is_sufficient(grid, stencil_width))
            {
                std::cout << "Error: Ghost layer not big enough. Either run with one_sided_BC=true or create a ghost "
                             "layer that has a width of least " << stencil_width << " grid node(s)" << std::endl;
                abort();
            }
        }
    }
    
    switch(order)
    {
        case 1:
        case 3:
        case 5:
            upwind_gradient<Field_in, Sign, Gradient_out>(grid, one_sided_BC, order);
            break;
        default:
            auto &v_cl = create_vcluster();
            if (v_cl.rank() == 0) std::cout << "Order of accuracy chosen not valid. Using default order 1." <<
                        std::endl;
            upwind_gradient<Field_in, Sign, Gradient_out>(grid, one_sided_BC, 1);
            break;
    }
    
    
}


#endif //OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_UPWIND_GRADIENT_HPP