wave_pressure_openmp.cc

#include <vtestbed/base/constants.hh>
#include <vtestbed/config/openmp.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/derivative.hh>
#include <vtestbed/core/grid_interpolation.hh>
#include <vtestbed/core/wave_pressure.hh>

using vtb::core::Policy;

namespace vtb {

    namespace core {

        template <class T>
        class Wave_pressure_solver_openmp: public Wave_pressure_solver<T> {

        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;

        protected:

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

            void
            calculate_pressure(
                const grid4& phi_grid,
                const array4& velocity_potential,
                const array3& wavy_surface,
                array4& pressure);

            void
            calculate_pressure_force(
                const grid4& pressure_grid,
                array4& pressure,
                panel_array& panels);

        };

    }

}

template <class T>
void
vtb::core::Wave_pressure_solver_openmp<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]);
    array4 phi(velocity_potential(_,rect[2],rect[0],rect[1]));
    array4 pressure(smaller_grid.shape());
    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_openmp<T>
::calculate_pressure(
    const grid4& grid,
    const array4& phi,
    const array3& zeta,
    array4& pressure
) {
    using blitz::pow2;
    using namespace vtb::core;
    using vtb::base::constants;
    const auto& delta = grid.delta();
    const auto g = this->g();
    constexpr auto p0 = constants<T>::atmospheric_pressure();
    parallel_for_loop<4>(
        grid.shape(),
        [&](const int4& idx) {
            const T z = grid(idx(1),1);
            const auto& grad = gradient(idx, delta, phi);
            pressure(idx) = p0-this->density() * (
                grad(0) + T(0.5)*(pow2(grad(1)) + pow2(grad(2)) + pow2(grad(3))) +
                g*z
            );
        }
    );
}

template <class T>
void
vtb::core::Wave_pressure_solver_openmp<T>
::calculate_pressure_force(
    const grid4& pressure_grid,
    array4& pressure,
    panel_array& panels
) {
    using blitz::Range;
    const auto _ = Range::all();
    const auto& new_pressure_grid = pressure_grid.select(1,2,3);
    Array<T,3> new_pressure(pressure(pressure_grid.num_patches(0),_,_,_));
    Linear_interpolation<T,3> interpolate{new_pressure, new_pressure_grid};
    const int npanels = panels.size();
    #if defined(VTB_WITH_OPENMP)
    #pragma omp parallel for
    #endif
    for (int i=0; i<npanels; ++i) {
        panel_type& panel = panels[i];
        const auto& centre = panel.centre();
        vec3 x{centre(2), centre(0), centre(1)};
        panel.force() = interpolate(x) * (-panel.normal()) * panel.area();
    }
}

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