lagrange_weights.hpp

Go to the documentation of this file.

#ifndef LAGRANGE_WEIGHTS_HPP
#define LAGRANGE_WEIGHTS_HPP

#define LAGRANGE_MAX_ORDER 8

#include "globals.hpp"

struct lagrange_weights{

    int idx = 0;             
    double w[LAGRANGE_MAX_ORDER] = { 0.0 }; 

    KOKKOS_INLINE_FUNCTION lagrange_weights(){}

    // Allow set in the constructor for legibility
    KOKKOS_INLINE_FUNCTION lagrange_weights(double var, double inv_d, int order){
        set(var, inv_d, order);
    }

    KOKKOS_INLINE_FUNCTION void set(double var, double inv_d, int order){
        double var_over_d, cell_loc, delta;

        var_over_d = (order > 0) ? var*inv_d/order : var*inv_d;
        idx = floor(var_over_d);
        cell_loc = var_over_d - idx; // Varies from 0.0 to 1.0
        if (order > 0){
            delta = 1.0/order;
        }
        else if(cell_loc > 0.5){ // For order 0 set idx to nearest grid point
            idx += 1;
        }

        // double w_sum = 0.0, m_sum = 0.0, e_sum = 0.0; // Temporary for debugging

        for (int kk = 0; kk <= order; kk++){
            w[kk] = 1.0;
            for (int jj = 0; jj <= order; jj++){
                if (kk != jj){
                    w[kk] *= (cell_loc - jj*delta)/(kk*delta - jj*delta);
                }
            }
            // w_sum += w[kk]; // Temporary for debugging
            // m_sum += w[kk]*(double)kk/order; // Temporary for debugging
            // e_sum += w[kk]*((double)kk/order)*((double)kk/order); // Temporary for debugging
        }

        // Temporary for debugging
        // double tol = 1.0e-10;
        // assert_XGC(1.0 - tol <= w_sum && w_sum <= 1.0 + tol,
        //            "\nERROR: lagrange weights density error\n");
        // assert_XGC(cell_loc - tol <= m_sum && m_sum <= cell_loc + tol,
        //            "\nERROR: lagrange weights momentum error\n");
        // assert_XGC(cell_loc*cell_loc - tol <= e_sum && e_sum <= cell_loc*cell_loc + tol,
        //            "\nERROR: lagrange weights energy error\n");
    }

    KOKKOS_INLINE_FUNCTION void set_idx(int new_idx){
        idx = new_idx;
    }

};

#endif