wave_pressure_opencl.cc

#include <vtestbed/base/constants.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/wave_pressure.hh>
#include <vtestbed/opencl/derivative.hh>
#include <vtestbed/opencl/opencl.hh>
#include <vtestbed/opencl/pipeline.hh>
#include <vtestbed/opencl/vector.hh>

using vtb::core::Policy;
using vtb::core::Wave_pressure_solver;
using vtb::opencl::Buffer;
using vtb::opencl::Context;
using vtb::opencl::Context_base;
using vtb::opencl::make_vector;

namespace vtb {

    namespace core {

        template <class T>
        class Wave_pressure_solver_opencl:
            public Wave_pressure_solver<T>,
            public Context_base {

        private:
            using base_type = Wave_pressure_solver<T>;

        public:
            using typename base_type::panel_type;
            using typename base_type::panel_array;
            using typename base_type::grid3;
            using typename base_type::grid4;
            using typename base_type::array3;
            using typename base_type::array4;
            using typename base_type::vec3;
            using derivative_type = vtb::opencl::Derivative<T,4>;

        private:
            derivative_type _derivative;
            Buffer<T> _d_derivative_phi[4];
            clx::kernel _calculate_pressure, _calculate_pressure_force;

        public:

            void
            calculate_pressure(
                const grid4& grid,
                const array4& velocity_potential,
                const array3& wavy_surface,
                Buffer<T>& pressure
            );

            void
            calculate_pressure_force(
                const grid4& pressure_grid,
                Buffer<T>& pressure,
                panel_array& panels
            );

            void solve(const grid3& wavy_surface_grid,
                       array3 wavy_surface,
                       const grid4& phi_grid,
                       array4 velocity_potential,
                       panel_array& panels) override;

            using Context_base::context;

            void context(Context* rhs) override {
                Context_base::context(rhs);
                auto prog = context()->compiler().compile("wave_pressure.cl");
                this->_calculate_pressure = prog.kernel("calculate_pressure");
                this->_calculate_pressure_force = prog.kernel("calculate_pressure_force");
                this->_derivative.context(rhs);
            }

        };

    }

}

template <class T>
void
vtb::core::Wave_pressure_solver_opencl<T>::solve(
    const grid3& wavy_surface_grid,
    array3 wavy_surface,
    const grid4& phi_grid,
    array4 velocity_potential,
    panel_array& panels
) {
    using blitz::Range;
    using blitz::RectDomain;
    grid3 grid_xyz = phi_grid.select(2,3,1);
    if (this->clip()) {
        grid_xyz = clamp(grid_xyz, panels);
        if (grid_xyz.empty()) { return; }
        grid_xyz.compact();
    }
    RectDomain<3> rect{grid_xyz.begin(), grid_xyz.end()-1};
    const auto& _ = Range::all();
    const grid4& smaller_grid = phi_grid(_,rect[2],rect[0],rect[1]);
    Array<T,4> phi(velocity_potential(_,rect[2],rect[0],rect[1]));
    Buffer<T> pressure;
    this->calculate_pressure(smaller_grid, phi, wavy_surface, pressure);
    this->calculate_pressure_force(smaller_grid, pressure, panels);
}

template <class T>
void
vtb::core::Wave_pressure_solver_opencl<T>
::calculate_pressure(
    const grid4& grid,
    const array4& phi,
    const array3& wavy_surface,
    Buffer<T>& d_pressure
) {
    using vtb::core::int1;
    using vtb::base::constants;
    const auto& shape = grid.shape();
    const size_t pressure_size = product(shape);
    auto context = this->context();
    auto& ppl = context->pipeline();
    Buffer<T> d_phi;
    ppl.copy(phi, d_phi);
    ppl.step();
    auto& d_derivative_phi = this->_d_derivative_phi;
    for (int i=0; i<4; ++i) {
        this->_derivative.solve(d_phi, grid, 0, d_derivative_phi[i]);
    }
    auto& kernel = this->_calculate_pressure;
    Buffer<T> d_wavy_surface;
    ppl.copy(wavy_surface, d_wavy_surface);
    ppl.allocate(pressure_size, d_pressure);
    kernel.arguments(
        d_derivative_phi[0], d_derivative_phi[1], d_derivative_phi[2], d_derivative_phi[3],
        d_wavy_surface,
        grid.lbound(1), grid.delta(1),
        shape(1), shape(2), shape(3),
        d_pressure);
    ppl.step();
    ppl.kernel(kernel, clx::range{pressure_size});
    ppl.step();
}

template <class T>
void
vtb::core::Wave_pressure_solver_opencl<T>::calculate_pressure_force(
    const grid4& pressure_grid,
    Buffer<T>& d_pressure,
    panel_array& panels
) {
    auto context = this->context();
    const auto& new_pressure_grid = pressure_grid.select(1,2,3);
    Buffer<panel_type> d_panels;
    auto& ppl = context->pipeline();
    ppl.copy(panels, d_panels);
    auto& kernel = this->_calculate_pressure_force;
    int offset = (pressure_grid.extent(0)-1)*product(new_pressure_grid.shape());
    kernel.argument(0, make_vector(new_pressure_grid.lbound()));
    kernel.argument(1, make_vector(new_pressure_grid.ubound()));
    kernel.argument(2, d_pressure);
    kernel.argument(3, offset);
    kernel.argument(4, make_vector(new_pressure_grid.shape()));
    kernel.argument(5, d_panels);
    ppl.step();
    ppl.kernel(kernel, clx::range(panels.size()));
    ppl.step();
    ppl.copy(d_panels, panels);
    ppl.wait();
}

template <>
auto
vtb::core::make_wave_pressure_solver<VTB_REAL_TYPE,Policy::OpenCL>()
-> std::unique_ptr<Wave_pressure_solver<VTB_REAL_TYPE>> {
    return std::unique_ptr<Wave_pressure_solver<VTB_REAL_TYPE>>(
        new Wave_pressure_solver_opencl<VTB_REAL_TYPE>
    );
}