gerstner_opencl.cc

#include <array>
#include <atomic>
#include <chrono>
#include <cstdio>

#include <vtestbed/base/constants.hh>
#include <vtestbed/base/for_loop.hh>
#include <vtestbed/config/openmp.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/bisection.hh>
#include <vtestbed/core/gerstner.hh>
#include <vtestbed/core/linear_wave.hh>
#include <vtestbed/core/ship_hull_panel.hh>
#include <vtestbed/core/types.hh>
#include <vtestbed/geometry/cubic_spline.hh>
#include <vtestbed/geometry/spline_surface.hh>
#include <vtestbed/geometry/tetrahedron.hh>
#include <vtestbed/opencl/opencl.hh>
#include <vtestbed/opencl/pipeline.hh>
#include <vtestbed/opencl/vector.hh>

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 Gerstner_solver_opencl:
            public Gerstner_solver<T, 3>,
            public Context_base {

        private:
            using base_type = Gerstner_solver<T, 3>;
            using typename base_type::vertex_field_3d;
            using typename base_type::scalar_field_3d;
            using typename base_type::panel_array;
            using typename base_type::panel_type;
            using typename base_type::vertex_type;
            using typename base_type::grid4;
            using typename base_type::grid3;
            using typename base_type::wave_type;
            using typename base_type::ship_type;

            struct small_panel {
                vertex_type centre{T{}}, normal{T{}};
                inline small_panel(const vertex_type& c, const vertex_type& n):
                centre(c), normal(n) {}
            };
            using small_panel_array = std::vector<small_panel>;

            static_assert(sizeof(small_panel) == sizeof(vertex_type)*2, "bad small_panel");
            static_assert(sizeof(wave_type) == 4*sizeof(T), "bad wave_type");

        private:
            clx::kernel _compute_forces;
            std::array<clx::kernel,(1<<5)> _generate_field;
            Buffer<wave_type> _d_waves;
            Buffer<panel_type> _d_panels;
            Buffer<small_panel> _d_small_panels;
            Buffer<vertex_type> _d_surface;

        public:

            using Context_base::context;

            void context(Context* rhs) override {
                Context_base::context(rhs);
                using clock = std::chrono::system_clock;
                using std::chrono::seconds;
                auto t0 = clock::now();
                size_t nkernels = this->_generate_field.size();
                std::atomic<int> count{};
                #if defined(VTB_WITH_OPENMP)
                #pragma omp parallel for schedule(dynamic,1)
                #endif
                for (size_t i=0; i<nkernels; ++i) {
                    bool diffraction = i & 1,
                         radiation = i & 2,
                         position = i & 4,
                         velocity_and_potential = i & 8,
                         finite_depth = i & 16;
                    auto cc = context()->compiler();
                    auto options = cc.options();
                    if (diffraction) { options += " -DVTB_DIFFRACTION"; }
                    if (radiation) { options += " -DVTB_RADIATION"; }
                    if (position) { options += " -DVTB_POSITION"; }
                    if (velocity_and_potential) {
                        options += " -DVTB_VELOCITY -DVTB_POTENTIAL";
                    }
                    if (finite_depth) { options += " -DVTB_FINITE_DEPTH"; }
                    else { options += " -DVTB_INFINITE_DEPTH"; }
                    cc.options(options);
                    auto prog = cc.compile("gerstner.cl");
                    this->_generate_field[i] = prog.kernel("generate_field");
                    auto cnt = ++count;
                    if (clock::now()-t0 > seconds(1)) {
                        std::fprintf(stderr, "%5d/%-5lu compile gerstner\n", cnt, nkernels);
                    }
                }
                auto& cc = context()->compiler();
                auto prog = cc.compile("gerstner.cl");
                this->_compute_forces = prog.kernel("compute_forces");
            }

            inline clx::kernel& generate_field_kernel(
                bool position,
                bool velocity_and_potential,
                bool finite_depth
            ) {
                size_t i = 0;
                if (this->diffraction()) { i |= 1; }
                if (this->radiation()) { i |= 2; }
                if (position) { i |= 4; }
                if (velocity_and_potential) { i |= 8; }
                if (finite_depth) { i |= 16; }
                return this->_generate_field[i];
            }

            void compute_forces(
                const ship_type& ship,
                const grid4& grid_tzxy,
                panel_array & wetted_panels
            ) override;

            void compute_positions(
                const ship_type& ship,
                const panel_array& panels,
                const grid3& grid_txy,
                vertex_field_3d& result
            ) override;

            void generate_field(
                const grid3& grid_zxy,
                T t,
                const panel_array& all_panels,
                const ship_type& ship,
                vertex_field_3d* position,
                vertex_field_3d* velocity=nullptr,
                scalar_field_3d* potential=nullptr
            );

        };

    }

}



template <class T>
void
vtb::core::Gerstner_solver_opencl<T>
::compute_forces(
    const ship_type& ship,
    const grid4 & grid_tzxy,
    panel_array & wetted_panels
) {
    // clamp grid to panels
    grid3 grid_zxy = grid_tzxy.select(1,2,3);
    if (this->clip()) { grid_zxy = clamp(grid_zxy, wetted_panels); }
    if (grid_zxy.ubound(2) > 0) { grid_zxy.ubound(2) = 0; }
    if (grid_zxy.empty()) { return; }
    grid_zxy.compact();
    this->_velocity_grid_zxy = grid_zxy;
    // compute velocity at grid points
    const T t = grid_tzxy.ubound(0);
    auto& velocity = this->_velocity;
    auto& potential = this->_potential;
    velocity.resize(grid_zxy.shape());
    potential.resize(grid_zxy.shape());
    generate_field(grid_zxy, t, wetted_panels, ship, nullptr, &velocity, &potential);
    auto& ppl = context()->pipeline();
    auto& kernel = this->_compute_forces;
    Buffer<panel_type> d_wetted_panels;
    Buffer<vertex_type> d_velocity;
    ppl.copy(wetted_panels, d_wetted_panels);
    ppl.copy(velocity, d_velocity);
    kernel.arguments(make_vector(grid_zxy.lbound()), make_vector(grid_zxy.ubound()),
                     make_vector(grid_zxy.shape()), d_velocity, d_wetted_panels);
    ppl.step();
    ppl.kernel(kernel, clx::range(wetted_panels.size()));
    ppl.step();
    ppl.copy(d_wetted_panels, wetted_panels);
    ppl.wait();
}

template <class T>
void
vtb::core::Gerstner_solver_opencl<T>
::compute_positions(
    const ship_type& ship,
    const panel_array& panels,
    const grid3 & grid_txy,
    vertex_field_3d& surface
) {
    const auto& _ = blitz::Range::all();
    const auto& grid_xy = grid_txy.select(1, 2);
    const T t = grid_txy.ubound(0);
    grid3 grid_zxy{Grid<T,1>{T{},T{},1}, grid_txy.select(1), grid_txy.select(2)};
    Array<vertex_type,3> position(grid_zxy.shape());
    generate_field(grid_zxy, t, panels, ship, &position, nullptr, nullptr);
    surface(grid_txy.end(0)-1,_,_) = position(0,_,_);
}

template <class T> void
vtb::core::Gerstner_solver_opencl<T>::generate_field(
    const grid3& grid_zxy,
    T t,
    const panel_array& all_panels,
    const ship_type& ship,
    vertex_field_3d* position,
    vertex_field_3d* velocity,
    scalar_field_3d* potential
) {
    const auto h = this->depth();
    const bool infinite_depth = is_positive_infinity(h);
    const bool waterline_only = this->waterline_only();
    auto& ppl = context()->pipeline();
    const auto& cc = context()->compiler();
    auto& waves = this->_waves;
    auto& d_waves = this->_d_waves;
    ppl.copy(waves, d_waves);
    auto& kernel = generate_field_kernel(position, velocity && potential, !infinite_depth);
    kernel.arguments(t, h, d_waves, int(waves.size()), make_vector(grid_zxy.lbound()),
                     make_vector(grid_zxy.delta()), make_vector(grid_zxy.shape()));
    size_t argument_index = 7;
    if (this->diffraction() || this->radiation()) {
        small_panel_array panels;
        panels.reserve(all_panels.size());
        for (const auto& panel : all_panels) {
            if ((waterline_only && panel.waterline()) || (!waterline_only && panel.wetted())) {
                panels.emplace_back(panel.centre(), panel.normal());
            }
        }
        auto& d_panels = this->_d_small_panels;
        int npanels = panels.size();
        if (npanels != 0) {
            ppl.copy(panels, d_panels);
        }
        kernel.argument(argument_index++, npanels);
        kernel.argument(argument_index++, d_panels);
        auto num_local_panels = std::numeric_limits<size_t>::max();
        const auto& device = cc.devices().front();
        auto wg = kernel.work_group(device);
        auto local_memory_size = device.local_memory_size();
        num_local_panels = std::max(size_t(1),
            std::min({wg.size, local_memory_size / sizeof(small_panel), panels.size()}));
        kernel.argument(argument_index++, int(num_local_panels));
        kernel.argument(argument_index++, clx::local<small_panel>(num_local_panels));
    }
    if (position && !velocity && !potential) {
        Buffer<vertex_type> d_result;
        ppl.allocate(position->shape(), d_result);
        kernel.argument(argument_index++, d_result);
        ppl.step();
        ppl.kernel(kernel, grid_zxy.shape());
        ppl.step();
        ppl.copy(d_result, *position);
        ppl.wait();
    } else if (!position && velocity && potential) {
        Array<Vector<T,4>,3> velocity_and_potential(velocity->shape());
        Buffer<Vector<T,4>> d_velocity_and_potential;
        ppl.allocate(velocity->shape(), d_velocity_and_potential);
        kernel.argument(argument_index++, d_velocity_and_potential);
        ppl.step();
        ppl.kernel(kernel, grid_zxy.shape());
        ppl.step();
        ppl.copy(d_velocity_and_potential, velocity_and_potential);
        ppl.wait();
        auto first = velocity->data();
        for (const auto& v : velocity_and_potential) {
            *first++ = vertex_type(v(0),v(1),v(2));
        }
        auto first2 = potential->data();
        for (const auto& v : velocity_and_potential) { *first2++ = v(3); }
    } else {
        throw std::invalid_argument("this combination of position, velocity, potential "
                                    "is not compiled");
    }
}

template <>
std::unique_ptr<vtb::core::Gerstner_solver<VTB_REAL_TYPE,3>>
vtb::core::make_gerstner_solver<VTB_REAL_TYPE,3,vtb::core::Policy::OpenCL>() {
    return std::unique_ptr<vtb::core::Gerstner_solver<VTB_REAL_TYPE,3>>(
        new Gerstner_solver_opencl<VTB_REAL_TYPE>
    );
}