velocity_potential_opencl.cc

#include <vtestbed/base/constants.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/math.hh>
#include <vtestbed/core/velocity_potential.hh>
#include <vtestbed/opencl/fourier_transform.hh>
#include <vtestbed/opencl/pipeline.hh>

using vtb::opencl::Buffer;
using vtb::opencl::Context;
using vtb::opencl::Context_base;

namespace vtb {

    namespace core {

        template <class T>
        class Linear_velocity_potential_solver_opencl:
            public Linear_velocity_potential_solver_base<T>,
            public Context_base
        {

        private:
            using C = std::complex<T>;
            using int2 = blitz::TinyVector<int,2>;
            using kernel_type = clx::kernel;
            using fft_type = vtb::opencl::Fourier_transform<C,2>;

        private:
            fft_type _fft;
            Array<T,3> _wf;
            Buffer<T> _d_window_function;
            Buffer<T> _d_phi;
            Grid<T,3> _wfgrid;
            clx::program _prog;

        public:

            using Context_base::context;

            void
            solve(
                const Grid<T,4>& grid_tzxy,
                const Array<T,3>& wavy_surface,
                Array<T,4>& velocity_potential
            ) override;

            void
            window_function(
                const Grid<T,2>& wngrid,
                const T z,
                Array2<T>& result
            );

            void context(Context* rhs) override {
                Context_base::context(rhs);
                this->_fft.context(rhs);
                fft_type::precompile({10,10}, rhs);
                auto prog = context()->compiler().compile("velocity_potential.cl");
                this->_prog = prog;
            }

        };

    }

}

typedef VTB_REAL_TYPE T;
using vtb::core::Policy;

namespace {

//  #define VTB_DEBUG_VELOCITY_POTENTIAL

    #if defined(VTB_DEBUG_VELOCITY_POTENTIAL)
    template <class T, int N>
    void
    dump(vtb::opencl::Buffer<T> d_x, const blitz::TinyVector<int,N>& shape,
         const char* name, int i=-1) {
        blitz::Array<T,N> x(shape);
        ppl.copy(d_x, x);
        ppl.wait();
        std::stringstream filename;
        filename << name;
        if (i >=0 ) {
            filename << i;
        }
        std::ofstream(filename.str()) << x;
    }
    #endif

    #if defined(VTB_DEBUG_VELOCITY_POTENTIAL)
    #define VTB_DUMP(x, shape) ::dump(x, shape, #x)
    #define VTB_DUMP2(x, shape, i) ::dump(x, shape, #x, i)
    #define VTB_DUMP3(x, shape, name, i) ::dump(x, shape, name, i)
    #else
    #define VTB_DUMP(x, shape)
    #define VTB_DUMP2(x, shape, i)
    #define VTB_DUMP3(x, shape, name, i)
    #endif

}

template <>
void
vtb::core::Linear_velocity_potential_solver_opencl<T>::solve(
    const Grid<T,4>& grid_tzxy,
    const Array<T,3>& zeta,
    Array<T,4>& phi
) {
    using blitz::Range;
    const int nt = grid_tzxy.extent(0);
    const Grid<T,2> grid_xy(grid_tzxy.select(2,3));
    const Range _ = Range::all();
    Array<C,3> sf(zeta.shape());
    this->second_function(grid_tzxy, zeta, sf);
    // update FFT wave table for the new shape
    {
        this->_fft.shape(grid_xy.shape());
        #if defined(VTB_DEBUG_VELOCITY_POTENTIAL)
        std::ofstream out{"fft.cl"};
        this->_fft.dump(out);
        #endif
    }
    auto wnrange = this->wave_number_grid();
    wnrange.end() = grid_xy.shape();
    Grid<T,3> new_wfgrid{grid_tzxy.select(1), wnrange.select(0), wnrange.select(1)};
    if (wnrange.empty() || near(wnrange.ubound(), T{0}, T{1e-3})) {
        phi = 0;
        return;
    }
    // limit wave number by the grid size
    Vector<T,2> wnmax = base::constants<T>::pi(4) / (grid_xy.ubound() - grid_xy.lbound());
    wnrange.ubound() = min(wnmax, wnrange.ubound());
    wnrange.ubound() = min(wnrange.lbound(), wnrange.ubound());
    auto& ppl = context()->pipeline();
    auto& d_window_function = this->_d_window_function;
    if (!new_wfgrid.near(this->_wfgrid, T{1e-3})) {
        this->_wfgrid = new_wfgrid;
        const auto& window_function_shape = this->_wfgrid.shape();
        this->_wf.resize(window_function_shape);
        ppl.allocate(this->_wfgrid.num_elements(), d_window_function);
        const T h = this->depth();
        bool infinite_depth = is_positive_infinity(h);
        auto lbound{this->_wfgrid.lbound()};
        auto delta{this->_wfgrid.delta()};
        auto kernel_name = infinite_depth ? "lvps_wf_infinite_depth" : "lvps_wf_finite_depth";
        auto kernel = this->_prog.kernel(kernel_name);
        kernel.argument(0, lbound(0));
        kernel.argument(1, delta(0));
        kernel.argument(2, lbound(1));
        kernel.argument(3, delta(1));
        kernel.argument(4, lbound(2));
        kernel.argument(5, delta(2));
        if (infinite_depth) {
            kernel.argument(6, d_window_function);
        } else {
            kernel.argument(6, h);
            kernel.argument(7, d_window_function);
        }
        ppl.kernel(kernel, window_function_shape);
        ppl.step();
    }
    auto& d_phi = this->_d_phi;
    ppl.allocate(phi.size(), d_phi);
    Buffer<C> d_fft_sf, d_result;
    ppl.allocate(this->_wfgrid.num_elements(), d_result);
    // for each time instant
    for (int i=0; i<nt; ++i) {
        Array<C,2> fft_sf(grid_xy.shape());
        fft_sf = sf(i,_,_);
        ppl.copy(fft_sf, d_fft_sf);
        ppl.step();
        this->_fft.forward(d_fft_sf);
        VTB_DUMP2(fft_sf, grid_xy.shape(), i);
        const auto& window_function_shape = this->_wfgrid.shape();
        {
            auto kernel = this->_prog.kernel("lvps_multiply");
            kernel.argument(0, d_fft_sf);
            kernel.argument(1, d_window_function);
            kernel.argument(2, d_result);
            ppl.step();
            ppl.kernel(kernel, window_function_shape);
        }
        ppl.step();
        this->_fft.backward(d_result, window_function_shape(0));
        {
            auto kernel = this->_prog.kernel("lvps_copy_back");
            kernel.argument(0, i);
            kernel.argument(1, d_result);
            kernel.argument(2, d_phi);
            ppl.step();
            ppl.kernel(kernel, window_function_shape);
        }
        ppl.wait();
    }
    ppl.copy(d_phi, phi);
    ppl.wait();
    VTB_DUMP3(phi, phi.shape(), "phi", -1);
}

template <>
auto
vtb::core::make_velocity_potential_solver<T,3,Policy::OpenCL>() ->
std::unique_ptr<Velocity_potential_solver<T,3>> {
    return std::unique_ptr<Velocity_potential_solver<T,3>>(
        new Linear_velocity_potential_solver_opencl<T>
    );
}

template class vtb::core::Linear_velocity_potential_solver_opencl<T>;