flux_surface_average.hpp

Go to the documentation of this file.

#ifndef FLUX_SURFACE_AVERAGE_HPP
#define FLUX_SURFACE_AVERAGE_HPP

#include "my_subview.hpp"

extern "C" void convert_grid_2_001d_wrap(double* input, double* output);
extern "C" void convert_001d_2_grid_wrap(double* input, double* output);

extern "C" int get_grid_npsi_surf();

// Apply flux-surface average on closed flux surfaces. Values on open surfaces are set to zero.
// The flux surface average function is currently in fortran,
// so the input must first be copied from device memory if using GPU
// Temporary arrays are passed as arguments to avoid reallocation if used repeatedly.
template<class Device>
void flux_surface_average(const View<double*, CLayout, Device>& input_view, const View<double*, CLayout, HostType>& output_view, const View<double*, CLayout, HostType>& tmp00_surf){
    Kokkos::deep_copy(output_view, input_view);
    convert_grid_2_001d_wrap(output_view.data(),tmp00_surf.data());
    convert_001d_2_grid_wrap(tmp00_surf.data(),output_view.data());
}

// Apply flux-surface average on closed flux surfaces. Values on open surfaces are set to zero.
// The flux surface average function is currently in fortran,
// so the input must first be copied from device memory if using GPU
// Temporary arrays are passed as arguments to avoid reallocation if used repeatedly.
template<class Device>
void flux_surface_average_in_place(const View<double*, CLayout, Device>& input_view, View<double*, CLayout, HostType>& tmp, View<double*, CLayout, HostType>& tmp00_surf){
    Kokkos::deep_copy(tmp, input_view);
    convert_grid_2_001d_wrap(tmp.data(),tmp00_surf.data());
    convert_001d_2_grid_wrap(tmp00_surf.data(),tmp.data());
    Kokkos::deep_copy(input_view, tmp);
}

// Apply flux-surface average without providing temporary arrays
template<class Device>
void flux_surface_average_in_place(const View<double*, CLayout, Device>& input_view){

    // Use these temporary host arrays
    View<double*, CLayout, HostType> tmp00_surf(NoInit("tmp00_surf"),get_grid_npsi_surf());
    View<double*, CLayout, HostType> tmp(NoInit("tmp"), input_view.layout());

    flux_surface_average_in_place(input_view, tmp, tmp00_surf);
}

// Apply flux-surface average on a collection of fields. Contiguous dimension of the input view
// must be the node index.
template<class Device>
void flux_surface_average_in_place(const View<double**, CLayout, Device>& input_view){
    // Infer problem size
    int n_arrays = input_view.extent(0);
    int n_nodes = input_view.extent(1);

    // Use these temporary host arrays
    View<double*, CLayout, HostType> tmp00_surf(NoInit("tmp00_surf"),get_grid_npsi_surf());
    View<double*, CLayout, HostType> tmp(NoInit("tmp"), n_nodes);

    for(int i=0; i<n_arrays; i++){
        auto one_field = my_subview(input_view,i);
        flux_surface_average_in_place(one_field, tmp, tmp00_surf);
    }
}

#endif