disk_spreading.H

Go to the documentation of this file.

#ifndef DISK_SPREADING_H_
#define DISK_SPREADING_H_
 
#include "AMReX_REAL.H"
#include "amr-wind/wind_energy/actuator/actuator_utils.H"
#include "amr-wind/wind_energy/actuator/disk/UniformCt.H"
#include "amr-wind/core/FieldRepo.H"
 
using namespace amrex::literals;
 
namespace amr_wind::actuator::ops {

template <typename T>
class SpreadingFunction
{
public:
    using OwnerType = T;
 
    void operator()(
        const T& actObj,
        const int lev,
        const amrex::MFIter& mfi,
        const amrex::Geometry& geom)
    {
        (this->*m_function)(actObj, lev, mfi, geom);
    }
 
    SpreadingFunction(const SpreadingFunction&) = delete;
    void operator=(const SpreadingFunction&) = delete;
 
    void (SpreadingFunction::*m_function)(
        const T& actObj,
        const int,
        const amrex::MFIter&,
        const amrex::Geometry&);
 
    void uniform_gaussian_spreading(
        const T& actObj,
        const int lev,
        const amrex::MFIter& mfi,
        const amrex::Geometry& geom)
    {
        const auto& bx = mfi.tilebox();
 
        const auto bxa = amr_wind::utils::realbox_to_box(
            actObj.m_data.info().bound_box, geom);
        const auto& bxi = bx & bxa;
        if (bxi.isEmpty()) {
            return;
        }
 
        const auto& sarr = actObj.m_act_src(lev).array(mfi);
        const auto& problo = geom.ProbLoArray();
        const auto& dx = geom.CellSizeArray();
 
        const auto& data = actObj.m_data.meta();
 
        const vs::Vector epsilon = vs::Vector::one() * data.epsilon;
        const vs::Vector m_normal(data.normal_vec);
        const auto* pos = actObj.m_pos.data();
        const auto* force = actObj.m_force.data();
        const int npts = data.num_force_pts;
        const int nForceTheta = data.num_force_theta_pts;
        const auto dTheta = ::amr_wind::utils::two_pi() / nForceTheta;
 
        amrex::ParallelFor(
            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                const vs::Vector cc{
                    problo[0] + (i + 0.5_rt) * dx[0],
                    problo[1] + (j + 0.5_rt) * dx[1],
                    problo[2] + (k + 0.5_rt) * dx[2],
                };
 
                amrex::RealArray src_force = {0.0_rt};
                for (int ip = 0; ip < npts; ++ip) {
                    const auto& pforce = force[ip] / nForceTheta;
                    const auto pLoc = pos[ip];
 
                    for (int it = 0; it < nForceTheta; ++it) {
                        const amrex::Real angle =
                            ::amr_wind::utils::degrees(it * dTheta);
                        const auto rotMatrix = vs::quaternion(m_normal, angle);
                        const auto diskPoint = pLoc & rotMatrix;
                        const auto distance = diskPoint - cc;
                        const auto projection_weight =
                            utils::gaussian3d(distance, epsilon);
 
                        src_force[0] += projection_weight * pforce.x();
                        src_force[1] += projection_weight * pforce.y();
                        src_force[2] += projection_weight * pforce.z();
                    }
                }
 
                sarr(i, j, k, 0) += src_force[0];
                sarr(i, j, k, 1) += src_force[1];
                sarr(i, j, k, 2) += src_force[2];
            });
    }
 
    void linear_basis_spreading(
        const T& actObj,
        const int lev,
        const amrex::MFIter& mfi,
        const amrex::Geometry& geom)
    {
        const auto& bx = mfi.tilebox();
 
        const auto bxa = amr_wind::utils::realbox_to_box(
            actObj.m_data.info().bound_box, geom);
        const auto& bxi = bx & bxa;
        if (bxi.isEmpty()) {
            return;
        }
 
        const auto& sarr = actObj.m_act_src(lev).array(mfi);
        const auto& problo = geom.ProbLoArray();
        const auto& dx = geom.CellSizeArray();
 
        const auto& data = actObj.m_data.meta();
 
        const amrex::Real dR = data.dr;
        const amrex::Real epsilon = data.epsilon;
        const vs::Vector m_normal(data.normal_vec);
        const vs::Vector m_origin(data.center);
        const auto* pos = actObj.m_pos.data();
        const auto* force = actObj.m_force.data();
        const int npts = data.num_force_pts;
 
        amrex::ParallelFor(
            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                const vs::Vector cc{
                    problo[0] + (i + 0.5_rt) * dx[0],
                    problo[1] + (j + 0.5_rt) * dx[1],
                    problo[2] + (k + 0.5_rt) * dx[2],
                };
 
                amrex::RealArray src_force = {0.0_rt};
                for (int ip = 0; ip < npts; ++ip) {
                    const auto R = utils::delta_pnts_cyl(
                                       m_origin, m_normal, m_origin, pos[ip])
                                       .x();
                    const auto dist_on_disk =
                        utils::delta_pnts_cyl(m_origin, m_normal, cc, pos[ip]);
                    const auto& pforce = force[ip];
 
                    const amrex::Real weight_R =
                        utils::linear_basis_1d(dist_on_disk.x(), dR);
                    const amrex::Real weight_T =
                        1.0_rt / (::amr_wind::utils::two_pi() * R);
                    const amrex::Real weight_N =
                        utils::gaussian1d(dist_on_disk.z(), epsilon);
                    const auto projection_weight =
                        weight_R * weight_T * weight_N;
 
                    src_force[0] += projection_weight * pforce.x();
                    src_force[1] += projection_weight * pforce.y();
                    src_force[2] += projection_weight * pforce.z();
                }
 
                sarr(i, j, k, 0) += src_force[0];
                sarr(i, j, k, 1) += src_force[1];
                sarr(i, j, k, 2) += src_force[2];
            });
    }
    void linear_basis_in_theta(
        const T& actObj,
        const int lev,
        const amrex::MFIter& mfi,
        const amrex::Geometry& geom)
    {
        const auto& bx = mfi.tilebox();
 
        const auto bxa = amr_wind::utils::realbox_to_box(
            actObj.m_data.info().bound_box, geom);
        const auto& bxi = bx & bxa;
        if (bxi.isEmpty()) {
            return;
        }
 
        const auto& sarr = actObj.m_act_src(lev).array(mfi);
        const auto& problo = geom.ProbLoArray();
        const auto& dx = geom.CellSizeArray();
 
        const auto& data = actObj.m_data.meta();
 
        const amrex::Real dR = data.dr;
        const amrex::Real dTheta =
            ::amr_wind::utils::two_pi() / data.num_vel_pts_t;
        const amrex::Real epsilon = data.epsilon;
        const vs::Vector m_normal(data.normal_vec);
        const vs::Vector m_origin(data.center);
        const auto* pos = actObj.m_pos.data();
        const auto* force = actObj.m_force.data();
        const int npts = data.num_force_pts;
 
        amrex::ParallelFor(
            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                const vs::Vector cc{
                    problo[0] + (i + 0.5_rt) * dx[0],
                    problo[1] + (j + 0.5_rt) * dx[1],
                    problo[2] + (k + 0.5_rt) * dx[2],
                };
 
                amrex::RealArray src_force = {0.0_rt};
                for (int ip = 0; ip < npts; ++ip) {
                    const auto radius =
                        utils::delta_pnts_cyl(
                            m_origin, m_normal, m_origin, pos[ip])
                            .x();
                    const auto dArc = radius * dTheta;
                    const auto dist_on_disk =
                        utils::delta_pnts_cyl(m_origin, m_normal, cc, pos[ip]);
                    const amrex::Real arclength = dist_on_disk.y() * radius;
                    const auto& pforce = force[ip];
 
                    const amrex::Real weight_R =
                        utils::linear_basis_1d(dist_on_disk.x(), dR);
                    const amrex::Real weight_T =
                        utils::linear_basis_1d(arclength, dArc);
                    const amrex::Real weight_N =
                        utils::gaussian1d(dist_on_disk.z(), epsilon);
                    const auto projection_weight =
                        weight_R * weight_T * weight_N;
 
                    src_force[0] += projection_weight * pforce.x();
                    src_force[1] += projection_weight * pforce.y();
                    src_force[2] += projection_weight * pforce.z();
                }
 
                sarr(i, j, k, 0) += src_force[0];
                sarr(i, j, k, 1) += src_force[1];
                sarr(i, j, k, 2) += src_force[2];
            });
    }
 
    SpreadingFunction() : m_function(&SpreadingFunction::linear_basis_spreading)
    {}
    void initialize(const std::string& key)
    {
        if (std::is_same<UniformCt, typename OwnerType::TraitType>::value) {
            if (key == "UniformGaussian") {
                m_function = &SpreadingFunction::uniform_gaussian_spreading;
            } else if (key == "LinearBasis") {
                m_function = &SpreadingFunction::linear_basis_spreading;
            } else {
                amrex::Abort("Invalid spreading type");
            }
        } else {
            m_function = &SpreadingFunction::linear_basis_in_theta;
        }
    }
};
} // namespace amr_wind::actuator::ops
#endif /* DISK_SPREADING_H_ */