wetted_surface_solver_opencl.cc

#include <vtestbed/config/opencl.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/policy.hh>
#include <vtestbed/core/ship.hh>
#include <vtestbed/core/ship_hull_panel.hh>
#include <vtestbed/core/wetted_surface_solver.hh>
#include <vtestbed/geometry/face.hh>
#include <vtestbed/opencl/opencl.hh>
#include <vtestbed/opencl/vector.hh>

//#define VTB_DEBUG_WETTED_PANELS_OPENCL

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, int N>
        class Wetted_surface_solver_opencl:
            public Wetted_surface_solver<T,N>,
            public Context_base {

        private:
            using base_type = Wetted_surface_solver<T,N>;
            using typename base_type::ship_type;
            using typename base_type::grid_type;
            using typename base_type::panel_type;
            using typename base_type::waterline_type;
            using typename base_type::vertex_type;
            //using vertex_type = vtb::geometry::Vertex<T,N>;
            using line_segment_type = typename waterline_type::value_type;
            using face_type = vtb::geometry::Face<N>;

        private:
            Buffer<vertex_type> _d_vertices;
            Buffer<face_type> _d_faces;
            Buffer<line_segment_type> _d_waterline;
            clx::kernel _kernel_regular, _kernel_irregular;

        public:

            using Context_base::context;

            void context(Context* rhs) override {
                Context_base::context(rhs);
                auto cc = context()->compiler();
                auto options = cc.options();
                cc.options(options + " -DVTB_REGULAR");
                this->_kernel_regular =
                    cc.compile("wetted_surface.cl").kernel("wetted_panels");
                cc.options(options + " -DVTB_IRREGULAR");
                this->_kernel_irregular =
                    cc.compile("wetted_surface.cl").kernel("wetted_panels");
            }

            void
            solve(
                const ship_type& ship,
                const grid_type& grid_txyz,
                Array3<T> wavy_surface
            ) override;

            void
            solve(
                const ship_type& ship,
                const grid_type& grid_txyz,
                Array3<vertex_type> wavy_surface
            ) override;

        private:

            void
            do_solve(const ship_type& ship,
                     const grid_type& grid_txyz,
                     const void* wavy_surface,
                     size_t wavy_surface_size,
                     clx::kernel& kernel);

        };

    }

}
template <class T, int N>
void
vtb::core::Wetted_surface_solver_opencl<T,N>::do_solve(
    const ship_type& ship,
    const grid_type& grid_txyz,
    const void* wavy_surface,
    size_t wavy_surface_size,
    clx::kernel& kernel
) {
    const auto& vertices = ship.hull().vertices();
    const auto& faces = ship.hull().faces();
    const int nvertices = vertices.size();
    const int nfaces = faces.size();
    auto& waterline = this->_waterline;
    auto& all_panels = this->_all_panels;
    auto& mask = this->_wetted_vertices_mask;
    mask.resize(nvertices);
    mask.assign(nvertices, false);
    all_panels.resize(nfaces);
    waterline.resize(nfaces);
    Grid<T,2> grid_xy = grid_txyz.select(1,2);
    #if defined(VTB_DEBUG_WETTED_PANELS_OPENCL)
    std::clog << "wavy surface size: " << grid_xy << std::endl;
    #endif
    int offset = (grid_txyz.extent(0)-1)*product(grid_xy.shape());
    Buffer<panel_type> d_panels;
    clx::buffer d_wavy_surface;
    auto& ppl = this->context()->pipeline();
    ppl.allocate(all_panels, d_panels);
    ppl.allocate(waterline, this->_d_waterline);
    ppl.copy(vertices, this->_d_vertices);
    ppl.copy(faces, this->_d_faces);
    ppl.copy(wavy_surface, wavy_surface_size, d_wavy_surface);
    kernel.argument(0, nfaces);
    kernel.argument(1, this->_d_faces);
    kernel.argument(2, this->_d_vertices);
    kernel.argument(3, d_panels);
    kernel.argument(4, make_vector(ship.position()));
    kernel.argument(5, make_vector(ship.rotation_matrix()));
    kernel.argument(6, d_wavy_surface);
    kernel.argument(7, offset);
    kernel.argument(8, make_vector(grid_xy.lbound()));
    kernel.argument(9, make_vector(grid_xy.ubound()));
    kernel.argument(10, make_vector(grid_xy.shape()));
    kernel.argument(11, this->_d_waterline);
    ppl.step();
    ppl.kernel(kernel, clx::range(nfaces));
    ppl.step();
    ppl.copy(d_panels, all_panels);
    ppl.copy(this->_d_waterline, waterline);
    ppl.wait();
    auto& dry_panels = this->_dry_panels;
    dry_panels.clear();
    dry_panels.reserve(nfaces/2);
    auto& wetted_panels = this->_wetted_panels;
    wetted_panels.clear();
    wetted_panels.reserve(nfaces/2);
    for (int i=0; i<nfaces; ++i) {
        const auto& panel = all_panels[i];
        if (panel.wetted()) {
            wetted_panels.emplace_back(panel);
            const auto& f = faces[i];
            mask[f[0]] = true, mask[f[1]] = true, mask[f[2]] = true;
        } else {
            dry_panels.emplace_back(panel);
        }
    }
    int new_waterline_size = 0;
    for (int i=0; i<nfaces; ++i) {
        const auto& segment = waterline[i];
        if (all(segment.front() == 0 && segment.back() == 0)) {
            new_waterline_size = i;
            break;
        }
    }
    waterline.resize(new_waterline_size);
}

template <class T, int N>
void
vtb::core::Wetted_surface_solver_opencl<T,N>::solve(
    const ship_type& ship,
    const grid_type& grid_txyz,
    Array3<vertex_type> wavy_surface
) {
    if (!wavy_surface.isStorageContiguous()) { wavy_surface.makeUnique(); }
    do_solve(ship, grid_txyz, wavy_surface.data(), wavy_surface.size()*sizeof(vertex_type),
             this->_kernel_irregular);
}

template <class T, int N>
void
vtb::core::Wetted_surface_solver_opencl<T,N>::solve(
    const ship_type& ship,
    const grid_type& grid_txyz,
    Array3<T> wavy_surface
) {
    if (!wavy_surface.isStorageContiguous()) { wavy_surface.makeUnique(); }
    do_solve(ship, grid_txyz, wavy_surface.data(), wavy_surface.size()*sizeof(T),
             this->_kernel_regular);
}

template <>
auto
vtb::core::make_wetted_surface_solver<VTB_REAL_TYPE,3,vtb::core::Policy::OpenCL>() ->
std::unique_ptr<Wetted_surface_solver<VTB_REAL_TYPE,3>> {
    return std::unique_ptr<Wetted_surface_solver<VTB_REAL_TYPE,3>>(
        new Wetted_surface_solver_opencl<VTB_REAL_TYPE,3>
    );
}

template class vtb::core::Wetted_surface_solver_opencl<VTB_REAL_TYPE,3>;