anlt_wind_solver_openmp.cc

#include <vtestbed/base/blitz.hh>
#include <vtestbed/base/constants.hh>
#include <vtestbed/config/openmp.hh>
#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/anlt_wind_solver.hh>
#include <vtestbed/core/policy.hh>
#include <vtestbed/core/ship_hull_panel.hh>
#include <vtestbed/geometry/tetrahedron.hh>
#include <vtestbed/linalg/linear_algebra.hh>

//#define VTB_DEBUG_WIND
//#define VTB_DEBUG_WIND_CYLINDER

#if defined(VTB_DEBUG_WIND)
#include <fstream>
#endif

namespace vtb {

    namespace core {

        template <class T>
        class ANLT_wind_solver_openmp: public ANLT_wind_solver<T> {

        private:
            using base_type = ANLT_wind_solver<T>;
            using typename base_type::array_type;
            using typename base_type::wind_field_type;
            using typename base_type::grid_type;
            using typename base_type::panel_type;
            using typename base_type::panel_array;
            using typename base_type::ray_type;
            using typename base_type::ray_array;

        public:
            void step(const grid_type& grid, panel_array& panels) override;

        };
    }
}

template <class T>
void vtb::core::ANLT_wind_solver_openmp<T>::step(
    const grid_type& grid,
    panel_array& panels
) {
    using vec3 = blitz::TinyVector<T,3>;
    using vtb::base::constants;
    using std::pow;
    vec3 u(0,50,0);
    auto& wind_field = this->_wind_field;
    wind_field.resize(grid.shape());
    wind_field = T{};
    auto stratification = this->stratification();
    const auto& delta = grid.delta();
    #if defined(VTB_DEBUG_WIND_CYLINDER)
    constexpr T R = T{54.8858}/2;
    T U = u(0);
    #endif
    if (this->near_the_boundary()) {
        parallel_for_loop<3>(
            grid.shape(),
            [&] (const int3& index) {
                const auto& r = grid(index, delta);
                vec3 sum{T{}};
                T total_weight{T{}};
                for (const auto& panel : panels) {
                    const auto& n = panel.normal();
                    const auto& S = panel.centre();
                    //if ((dot(n,u) > 0) || (dot(r-S, n) < 0)) {
                    //  continue;
                    //}
                    //if (dot(n,u) > 0) { continue; }
                    if (stratification) {
                        T coef = 25.0;
                        T alpha = 1.15;
                        T z0 = 10.0;
                        T S_z = panel.centre()(2);
                        T z = r(2);
                        vec3 u_r = u - 2 * dot(u,n) * n;
                        vec3 u1(0, 0, (1+alpha)*pow(S_z,alpha)*dot(u, S));
                        vec3 u2(0, 0, (1+alpha)*pow(S_z,alpha));
                        vec3 u3(1, 1, (1+alpha)*pow(S_z,alpha + 1));
                        vec3 u4(0, 0, (1+alpha)*pow(z,alpha));
                        T S_r = sqrt(dot(S-r,S-r));
                        vec3 delta_phi(0, 0, 0);
                        T z_a = pow(z0, -alpha)/(alpha + 1);

                        //T C2 = 1.0;
                        T C2 = dot(vec3(u*pow(S_z,alpha+1)+ u1), n)/
                            dot(vec3(u2*dot(u_r,S) + u3*u_r), n);
                        //delta_phi += z_a*u;
                        //delta_phi = z_a*C2*u_r*exp(-S_r/coef)*pow(z,1+alpha)+
                        delta_phi += z_a*C2*u_r*exp(-S_r/coef)*pow(z,1+alpha);
                        //delta_phi += z_a*C2*u_r*exp(-S_r/coef)*pow(z,1+alpha);
                            //z_a*C2*dot(u_r,r)*exp(-S_r);
                            //z_a*C2*dot(u_r,r)*exp(-S_r)*
                            //(u2 + (S-r)*pow(z,1+alpha))/S_r;

                        sum += delta_phi;
                    } else {
                        vec3 u_r = u - T{2}*dot(u,n)*n;
                        T weight = T{1} / (T{1}+sum_of_squares(vec3(r-S)));
                        sum += u_r * weight;
                        total_weight += weight;
                    }
                }
                if (!all(isfinite(sum))) {
                    throw std::runtime_error("wind field is not finite");
                }
                if (total_weight != 0) { sum /= total_weight; }
                if (stratification) {
                    T alpha = 1.15;
                    T z_a = pow(10.0, -alpha)/(alpha + 1);
                    //T z_a = pow(10.0, -alpha);
                    sum += u * z_a * pow(r(2), alpha);
                    //sum += u * z_a;
                } else {
                    sum += u;
                }
                #if defined(VTB_DEBUG_WIND_CYLINDER)
                // solution for potential flow around a cylinder
                auto x = r(0), y = r(1);
                auto pow2 = [] (T x) { return x*x; };
                sum = vec3{U*(R*R*(y*y-x*x)+pow2(x*x+y*y)), -2*R*R*U*x*y, 0}/pow2(x*x+y*y);
                #endif
                wind_field(index) = sum;
            }
        );
    }
    #if defined(VTB_DEBUG_WIND)
    {
        #if defined(VTB_DEBUG_WIND_CYLINDER)
        std::ofstream out("wind-orig");
        #else
        std::ofstream out("wind-our");
        #endif
        for_loop<3>(
            grid.shape(),
            [&] (const int3& index) {
                if (index(2) != grid.extent(2)/2) { return; }
                const auto& v = wind_field(index);
                const auto& r = grid(index, delta);
                #if defined(VTB_DEBUG_WIND_CYLINDER)
                if (r(0)*r(0)+r(1)*r(1) < R*R) { return; }
                #endif
                out << r(0) << ' ' << r(1) << ' ';
                out << v(0) << ' ' << v(1) << ' ';
                out << '\n';
            }
        );
    }
    auto m = minmax(wind_field);
    std::clog << "phi_min=" << m.min << std::endl;
    std::clog << "phi_max=" << m.max << std::endl;
    #endif
    auto& wind_field_boundary = this->_wind_field_boundary;
    if (this->on_the_boundary()) {
        const int npanels = panels.size();
        wind_field_boundary.resize(npanels);
        constexpr auto rho = constants<T>::air_density();
        constexpr auto g = constants<T>::g();
        constexpr auto p0 = constants<T>::atmospheric_pressure();
        #if defined(VTB_WITH_OPENMP)
        #pragma omp parallel for
        #endif
        for (int i=0; i<npanels; ++i) {
            auto& panel = panels[i];
            const auto& r = panel.centre();
            const auto& n = panel.normal();
            vec3 u_r = u - 2*dot(u,n)*n;
            vec3 velocity = u + u_r;
            auto z = r(2);
            auto pressure = p0-rho*(T{0.5}*sum_of_squares(velocity) + g*z);
            panel.force() = pressure * (-n) * panel.area();
            wind_field_boundary[i] = ray_type(r, velocity);
        }
    } else {
        wind_field_boundary.clear();
    }
}

template <>
auto
vtb::core::make_anlt_wind_solver<VTB_REAL_TYPE,vtb::core::Policy::OpenMP>() ->
std::unique_ptr<ANLT_wind_solver<VTB_REAL_TYPE>> {
    return std::unique_ptr<ANLT_wind_solver<VTB_REAL_TYPE>>(
        new ANLT_wind_solver_openmp<VTB_REAL_TYPE>
    );
}

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