core/testbed.cc

#include <cstring>
#include <fstream>
#include <iostream>
#include <limits>

#include <vtestbed/config/opencl.hh>
#include <vtestbed/config/openmp.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/anlt_wind_solver.hh>
#include <vtestbed/core/arma_wind_solver.hh>
#include <vtestbed/core/cosine_wave.hh>
#include <vtestbed/core/propulsor.hh>
#include <vtestbed/core/ship_motion_solver.hh>
#include <vtestbed/core/stokes_wave.hh>
#include <vtestbed/core/testbed.hh>
#include <vtestbed/core/wetted_surface_solver.hh>
#include <vtestbed/geometry/basis.hh>
#include <vtestbed/profile/profile.hh>

#if defined(VTB_WITH_OPENCL)
#include <vtestbed/opencl/opencl.hh>
#endif

static_assert(
    sizeof(vtb::core::Ship_hull_panel<VTB_REAL_TYPE,3>) ==
    sizeof(vtb::core::Vector3<VTB_REAL_TYPE>)*4,
    "bad size"
);

template <class T, int N>
vtb::core::Testbed<T,N>::Testbed():
_grid_txyz{
    {T{-1}, T(-100), T(-100), T{-100}},
    {T{0}, T(100), T(100), T{2}},
    {3, 128, 128, 64} // we need three points for zeta time derivative
},
_wind_grid{
    {-100/2,-100/2,0},
    {100/2,100/2,25},
    {40,40,10}
} {
    this->_wpdomain_xyz.lbound() = 0;
    this->_wpdomain_xyz.ubound() = 0;
    grid3 wsgrid{this->wavy_surface_grid()};
    this->_wavysurface.resize(wsgrid.shape());
    this->_wavysurface = 0;
    this->_irregular_wavysurface.resize(wsgrid.shape());
    parallel_for_loop<N>(
        wsgrid.shape(),
        [&] (const blitz::TinyVector<int,N>& idx) {
            auto txy = wsgrid(idx);
            this->_irregular_wavysurface(idx) = {txy(1), txy(2), 0};
        }
    );
    grid4 vpgrid{this->_grid_txyz.select(0,3,1,2)};
    this->_vpotential.resize(vpgrid.shape());
    this->_statistics_log.open("statistics.log");
    this->_statistics.write_header(this->_statistics_log);
    //this->_arma_wind_solver.reset(new arma_wind_solver_type);
    //this->_arma_wind_solver->initialize();
    this->_ship_motion_solver.reset(new ship_motion_solver_type);
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::step(T dt) {
    VTB_PROFILE_BLOCK("step");
    const T t0 = this->_grid_txyz.ubound(0);
    this->_grid_txyz.lbound(0) = t0;
    this->_grid_txyz.ubound(0) = t0 + dt;
    switch (this->get_solver_type()) {
        case Solver_type::regular_solver: this->step_regular(); break;
        case Solver_type::gerstner_solver: this->step_irregular(); break;
    }
    /*
    std::ofstream out("waterline");
    waterline().gnuplot(out);
    */
    //TODO:
    //this->_arma_wind_solver->step(this->_wssolver->dry_panels());
    //this->_wind_field.reference(this->_arma_wind_solver->wind_field());
    this->compute_wind();
    this->move_ship();
    this->record_statistics(this->_statistics);
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::step_regular() {
    this->generate_waves();
    this->compute_wetted_panels_regular();
    if (has_wetted_panels()) {
        this->clamp_grid_to_wetted_panels();
        if (has_wave_pressure_domain()) {
            this->compute_wave_numbers();
            this->compute_velocity_potential();
            this->compute_wave_pressure();
        }
    }
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::step_irregular() {
    {
        VTB_PROFILE_BLOCK("fluid_motion");
        this->_gerstner_solver->compute_positions(this->_ship, this->_wssolver->all_panels(),
            this->_grid_txyz.select(0,1,2), this->_irregular_wavysurface);
    }
    this->compute_wetted_panels_irregular();
    if (has_wetted_panels()) {
        this->clamp_grid_to_wetted_panels();
        if (has_wave_pressure_domain()) {
            VTB_PROFILE_BLOCK("compute_forces");
            // TODO implement for irregular wavy surface
            //this->compute_wave_numbers();
            this->_gerstner_solver->compute_forces(
                this->_ship,
                this->_grid_txyz.select(0,3,1,2),
                this->_wssolver->wetted_panels());
            this->_water_velocity.reference(this->_gerstner_solver->velocity());
        }
    }
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::reset() {
    this->_grid_txyz.lbound(0) = -1;
    this->_grid_txyz.ubound(0) = 0;
    this->_ship.reset();
    this->_statistics.clear();
    this->_statistics_log.close();
    this->_statistics_log.open("statistics.log");
    this->_wavysurface = 0;
    this->_vpotential = 0;
    this->_wind_field = T{};
    this->_irregular_wavysurface = T{};
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::generate_waves() {
    VTB_PROFILE_BLOCK("waves");
    grid3 wsgrid{this->wavy_surface_grid()};
    this->_wavysurface.reference(this->_wsgenerator->generate(wsgrid));
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::compute_wind() {
    VTB_PROFILE_BLOCK("wind");
    this->_anlt_wind_solver->step(this->_wind_grid, this->_wssolver->dry_panels());
    this->_wind_field.reference(this->_anlt_wind_solver->wind_field());
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::compute_wetted_panels_regular() {
    VTB_PROFILE_BLOCK("wetted_panels");
    this->_wssolver->solve(ship(), grid_txyz(), this->wavy_surface());
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::compute_wetted_panels_irregular() {
    VTB_PROFILE_BLOCK("wetted_panels");
    this->_wssolver->solve(ship(), grid_txyz(), this->irregular_wavy_surface());
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::clamp_grid_to_wetted_panels() {
    using vtb::base::MinMax;
    const auto& wetted_panels = this->_wssolver->wetted_panels();
    MinMax<T> mm[3];
    for (const auto& panel : wetted_panels) {
        const auto& centre = panel.centre();
        mm[0].update(centre(0));
        mm[1].update(centre(1));
        mm[2].update(centre(2));
    }
    vec3 x_min{mm[0].min(),mm[1].min(),mm[2].min()};
    vec3 x_max{mm[0].max(),mm[1].max(),mm[2].max()};
    // set Z dimension to the wetted panels Z range
    if (!wetted_panels.empty()) {
        T z_max = std::min(x_max(2), blitz::max(this->_wavysurface)) + 1;
        //T z_max = std::min(x_max(2), blitz::max(this->_wavysurface)) + 1;
        this->_grid_txyz.ubound(3) = z_max;
        this->_grid_txyz.lbound(3) = std::min(x_min(2)-1, z_max);
    }
    grid3 grid_xyz{this->_grid_txyz.select(1,2,3)};
    x_min = clamp(x_min, grid_xyz);
    x_max = clamp(x_max, grid_xyz);
    int3 idx_min = int3{blitz::floor(grid_xyz.index(x_min))};
    int3 idx_max = int3{blitz::ceil(grid_xyz.index(x_max))};
    idx_min = blitz::max(0, idx_min);
    idx_max = blitz::min(grid_xyz.shape()-1, idx_max);
    // Make X and Y extent a power of 2. Do not touch Z dimension.
    int3 extent = idx_max - idx_min + 1;
    if (blitz::all(extent > 0)) {
        for (int i=0; i<2; ++i) {
            int e = extent(i);
            if (!blitz::is_power_of_two(e)) {
                extent(i) = blitz::next_power_of_two(e);
            }
        }
        idx_max = idx_min + extent - 1;
        for (int i=0; i<2; ++i) {
            int diff = idx_max(i) - (grid_xyz.extent(i) - 1);
            if (diff > 0) {
                idx_min(i) -= diff;
                idx_max(i) -= diff;
            }
        }
        if (blitz::any(idx_min < 0)) {
            std::clog
                << "Please, increase wavy surface size to speed-up Fourier transforms"
                << std::endl;
        }
        extent = idx_max - idx_min + 1;
    }
    domain3 d{idx_min, idx_max};
    this->_wpdomain_xyz = d;
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::compute_wave_numbers() {
    const auto& _ = blitz::Range::all();
    const auto& grid_xy = this->_grid_txyz.select(1,2);
    const auto& zeta = this->_wavysurface;
    const int slice_t = zeta.extent(0)-1;
    const auto& w = wave_statistics(zeta(slice_t,_,_), grid_xy);
    auto& waves = this->_statistics.waves();
    waves(1) = w(0);
    waves(2) = w(1);
}

template <class T, int N>
auto
vtb::core::Testbed<T,N>::wave_number_grid_2d() const -> grid2 {
    return wave_number_grid(this->_statistics.waves()).select(1,2);
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::compute_velocity_potential() {
    VTB_PROFILE_BLOCK("wave_velocity");
    const auto _ = blitz::Range::all();
    _vpgrid_tzxy = _grid_txyz.select(0,3,1,2);
    _vpgrid_tzxy = _vpgrid_tzxy(_,_wpdomain_xyz[2],_wpdomain_xyz[0],_wpdomain_xyz[1]);
    _vpotential.resize(_vpgrid_tzxy.shape());
    _vpsolver->wave_number_grid(this->wave_number_grid_2d());
    _vpsolver->solve(
        _vpgrid_tzxy,
        _wavysurface(_,_wpdomain_xyz[0],_wpdomain_xyz[1]),
        _vpotential
    );
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::compute_wave_pressure() {
    VTB_PROFILE_BLOCK("wave_pressure");
    grid3 wsgrid{_vpgrid_tzxy.select(0,2,3)};
    this->_wpsolver->solve(
        wsgrid,
        this->_wavysurface,
        _vpgrid_tzxy,
        this->_vpotential,
        this->_wssolver->wetted_panels()
    );
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::move_ship() {
    VTB_PROFILE_BLOCK("ship_motions");
    const T t0 = this->_grid_txyz.lbound(0);
    const T t1 = this->_grid_txyz.ubound(0);
    this->_ship_motion_solver->solve(
        this->_ship,
        this->wetted_panels(),
        this->dry_panels(),
        this->waterline(),
        t0,
        t1
    );
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::record_statistics(statistics_type& stats) {
    using Record = typename statistics_type::Record;
    VTB_PROFILE_BLOCK("statistics");
    T t = this->_grid_txyz.ubound(0);
    const auto& ship = this->_ship;
    stats.record(Record::Time, t);
    stats.record(Record::Euler_angles, ship.euler_angles());
    stats.record(Record::Position, ship.position());
    stats.record(Record::Velocity, ship.velocity());
    stats.record(Record::Angular_velocity, vector_to(ship.rotation_matrix(), ship.angular_velocity()));
    stats.record(Record::Acceleration, vector_to(ship.rotation_matrix(), ship.acceleration()));
    stats.record(Record::Angular_acceleration, ship.angular_acceleration());
    stats.record(Record::Underwater_volume,
                 this->_ship_motion_solver->equation().displacement());
    stats.record(Record::Angular_momentum, ship.angular_momentum());
    const auto& zeta = this->_wavysurface;
    const auto nx = zeta.extent(1);
    const auto ny = zeta.extent(2);
    stats.record(Record::Wavy_surface, zeta(zeta.extent(0)-1,nx/2,ny/2));
    stats.record_waves();
    stats.write_back(this->_statistics_log);
    this->_statistics_log.flush();
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::openmp_solvers() {
    #if defined(VTB_WITH_OPENMP)
    auto solver_t = get_solver_type();
    this->velocity_potential_solver(make_velocity_potential_solver<T,3,Policy::OpenMP>());
    this->wave_pressure_solver(make_wave_pressure_solver<T,Policy::OpenMP>());
    this->wetted_surface_solver(make_wetted_surface_solver<T,3,Policy::OpenMP>());
    this->anlt_wind_solver(make_anlt_wind_solver<T,Policy::OpenMP>());
    using default_wave = Propagating_stokes_wave<T,3>;
    this->wavy_surface_generator(
        make_plain_wave_generator<T,3,default_wave,Policy::OpenMP>()
    );
    this->gerstner_solver(make_gerstner_solver<T,3,Policy::OpenMP>());
    set_solver_type(solver_t);
    #endif
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::opencl_solvers() {
    #if defined(VTB_WITH_OPENCL)
    auto solver_t = get_solver_type();
    this->velocity_potential_solver(make_velocity_potential_solver<T,3,Policy::OpenCL>());
    this->wave_pressure_solver(make_wave_pressure_solver<T,Policy::OpenCL>());
    this->wetted_surface_solver(make_wetted_surface_solver<T,3,Policy::OpenCL>());
    this->anlt_wind_solver(make_anlt_wind_solver<T,Policy::OpenCL>());
    using default_wave = Propagating_stokes_wave<T,3>;
    this->wavy_surface_generator(
        make_plain_wave_generator<T,3,default_wave,Policy::OpenCL>()
    );
    this->gerstner_solver(make_gerstner_solver<T,3,Policy::OpenCL>());
    set_solver_type(solver_t);
    #endif
}

template <class T, int N>
void
vtb::core::Testbed<T,N>::opencl_context(void* ptr) {
    #if defined(VTB_WITH_OPENCL)
    using namespace vtb::opencl;
    auto context = reinterpret_cast<Context*>(ptr);
    if (auto vpsolver = dynamic_cast<Context_base*>(this->_vpsolver.get())) {
        vpsolver->context(context);
    }
    if (auto wpsolver = dynamic_cast<Context_base*>(this->_wpsolver.get())) {
        wpsolver->context(context);
    }
    if (auto wssolver = dynamic_cast<Context_base*>(this->_wssolver.get())) {
        wssolver->context(context);
    }
    if (auto wsgenerator = dynamic_cast<Context_base*>(this->_wsgenerator.get())) {
        wsgenerator->context(context);
    }
    if (auto gerstner_solver = dynamic_cast<Context_base*>(this->_gerstner_solver.get())) {
        gerstner_solver->context(context);
    }
    if (auto wind_solver = dynamic_cast<Context_base*>(this->_anlt_wind_solver.get())) {
        wind_solver->context(context);
    }
    #endif
}

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