gerstner_openmp.cc

#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/gerstner.hh>
#include <vtestbed/core/grid_interpolation.hh>
#include <vtestbed/core/ship_hull_panel.hh>
#include <vtestbed/core/types.hh>
#include <vtestbed/geometry/linear_interpolation.hh>

namespace vtb {

    namespace core {

        template <class T>
        class Gerstner_solver_openmp: public Gerstner_solver<T,3> {

        protected:
            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::C;
            using typename base_type::ship_type;

            enum class Quantity { Position, Velocity };

        public:

            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_xyz,
                vertex_field_3d& result,
                T t,
                const panel_array& panels,
                const ship_type& ship,
                Quantity q,
                scalar_field_3d* potential=nullptr
            );

        };

    }

}

template <class T> void
vtb::core::Gerstner_solver_openmp<T>::compute_forces(
    const ship_type& ship,
    const grid4& grid_tzxy,
    panel_array & wetted_panels
) {
    int npanels = wetted_panels.size();
    // 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());
    auto old_diffraction = this->diffraction();
    this->diffraction(false);
    generate_field(grid_zxy, velocity, t, wetted_panels, ship,
                            Quantity::Velocity, &potential);
    this->diffraction(old_diffraction);
    // interpolate velocity and compute wave pressure
    Linear_function_irregular<T,vertex_type> function(velocity, grid_zxy);
    Grid_interpolation<T,3,Linear_function_irregular<T,vertex_type>,vertex_type>
        interpolate(grid_zxy, function);
    constexpr const auto g = vtb::base::constants<T>::g();
    constexpr const auto rho = vtb::base::constants<T>::water_density();
    constexpr const auto p0 = vtb::base::constants<T>::atmospheric_pressure();
    #if defined(VTB_WITH_OPENMP)
    #pragma omp parallel for
    #endif
    for (int i=0; i<npanels; ++i) {
        auto& panel = wetted_panels[i];
        const auto& centre = panel.centre();
        auto z = centre(2);
        const auto& velocity = interpolate({centre(2),centre(0),centre(1)});
        T pressure = p0 - rho*(T{0.5}*sum_of_squares(velocity) + g*z);
        panel.force() = pressure * (-panel.normal()) * panel.area();
    }
}

template <class T> void
vtb::core::Gerstner_solver_openmp<T>::generate_field(
    const grid3& grid_zxy,
    vertex_field_3d& result,
    T t,
    const panel_array& all_panels,
    const ship_type& ship,
    Quantity q,
    scalar_field_3d* potential
) {
    using std::cosh;
    using std::exp;
    using std::sqrt;
    using vtb::geometry::unit;
    using vec3 = Vector<T,3>;
    constexpr const auto pi = vtb::base::constants<T>::pi();
    constexpr const auto g = vtb::base::constants<T>::g();
    constexpr const C i(0,1);
    const auto h = this->depth();
    const bool infinite_depth = is_positive_infinity(h);
    const bool diffraction = this->diffraction();
    const bool radiation = this->radiation();
    const bool waterline_only = this->waterline_only();
    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);
        }
    }
    parallel_for_loop<3>(
        grid_zxy.shape(),
        [&] (const int3& idx) {
            const auto& zxy = grid_zxy(idx);
            const auto x = zxy(1), y = zxy(2), z = zxy(0);
            vertex_type sum_position{x,y,z};
            vertex_type sum_velocity{0,0,0};
            T sum_potential{};
            for (const auto& wave : this->_waves) {
                const auto& _k_ = wave.wave_number();
                T u = _k_(0), v = _k_(1);
                T k = sqrt(u*u+v*v);
                T omega = wave.angular_frequency();
                T ampl = wave.amplitude();
                T exp1 = infinite_depth ? exp(k*z) : cosh(k*(z+h));
                C exp2 = exp(i*(u*x + v*y - omega*t));
                C w = exp1 * exp2;
                vec3 sum_dr{T{}};
                if (diffraction) {
                    C exp3{};
                    vec3 di(u,v,0), zeta(x,y,z);
                    T sum_weight = 0;
                    for (const auto& panel : panels) {
                        const auto& centre = panel.centre();
                        const auto& n = panel.normal();
                        vec3 ds = 2*dot(di,n)*n;
                        vec3 dr = di-ds;
                        T weight = 1/(1+distance(zeta,centre));
                        //T weight = 1;
                        exp3 += exp(i*(dot(dr,zeta) + dot(ds,centre)))*weight;
                        //T{0.5}*(1+dot(unit(di),n));
                        sum_dr += dr;
                        sum_weight += weight;
                    }
                    if (sum_weight != 0) {
                        exp3 /= sum_weight, sum_dr /= sum_weight;
                    }
                    w += exp1*exp3*exp(-i*omega*t);
                }
                vec3 sum_wave_vector{T{}};
                if (radiation) {
                    C exp4{};
                    vec3 zeta(x,y,z);
                    T sum_weight = 0;
                    for (const auto& panel : panels) {
                        const auto& n = panel.normal();
                        const auto& centre = panel.centre();
                        vec3 velocity = ship.velocity() +
                            cross(ship.angular_velocity(), vec3(centre-ship.position()));
                        auto old_velocity = velocity;
                        velocity = dot(velocity,n)*n;
                        //if (dot(vec3(zeta-ship.position()),n) < 0) { velocity -= n; }
                        T weight = 1/(1+length(vec3(zeta-centre)));
                        //vec3 wave_vector = length(velocity) / g * n;
                        vec3 wave_vector = g / (length(velocity)) * n/2;
                        //vec3 wave_vector = dot(old_velocity,n)/2*n;
                        T omega = dot(wave_vector,velocity);
                        exp4 += exp(i*(dot(wave_vector,zeta) - omega*t))*weight;
                        sum_wave_vector += wave_vector;
                        sum_weight += weight;
                    }
                    if (sum_weight != 0) {
                        exp4 /= sum_weight, sum_wave_vector /= sum_weight;
                    }
                    w += exp1*exp4;
                }
                sum_potential += real(w);
                T scale = ampl*2*pi;
                if (q == Quantity::Position) {
                    sum_position(0) += (u+sum_dr(0)+sum_wave_vector(0))*real(i*w);
                    sum_position(1) += (v+sum_dr(1)+sum_wave_vector(1))*real(i*w);
                    sum_position(2) += scale*real((k+i*sum_dr(2)+i*sum_wave_vector(2))*w);
                } else {
                    sum_velocity(0) += omega*(u+sum_dr(0)+sum_wave_vector(0))*real(w);
                    sum_velocity(1) += omega*(v+sum_dr(1)+sum_wave_vector(1))*real(w);
                    sum_velocity(2) += -scale*omega*k*real(i*w);
                }
            }
            result(idx) = (q == Quantity::Position) ? sum_position : sum_velocity;
            if (potential) { (*potential)(idx) = sum_potential; }
        }
    );
}

template <class T> void
vtb::core::Gerstner_solver_openmp<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{
        {T{},grid_txy.lbound(1),grid_txy.lbound(2)},
        {T{},grid_txy.ubound(1),grid_txy.ubound(2)},
        {1,grid_txy.extent(1),grid_txy.extent(2)}};
    Array<vertex_type,3> position(grid_zxy.shape());
    generate_field(grid_zxy, position, t, panels, ship, Quantity::Position);
    surface(grid_txy.end(0)-1,_,_) = position(0,_,_);
}

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

template class vtb::core::Gerstner_solver_openmp<VTB_REAL_TYPE>;