Program Listing for File biot_savart.h

↰ Return to documentation for file (source/adem/biot_savart.h)

/*
 * biot_savart.h A parallelised biot-savart function for use with the Attached-Detached Eddy Method
 *
 * Author:              Tom Clark  (thclark @ github)
 *
 * Copyright (c) 2014-2019 Octue Ltd. All Rights Reserved.
 *
 */

#ifndef ES_FLOW_BIOT_SAVART_H
#define ES_FLOW_BIOT_SAVART_H

#include <tbb/tbb.h>
#include <tbb/blocked_range.h>
#include <tbb/parallel_for.h>
#include <tbb/atomic.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <iostream>
#include <Eigen/Dense>


using namespace Eigen;
using namespace std;


#define pi          3.14159265358979323846264;
#define fourpi      12.5663706143592172953851;
#define invfourpi   0.079577471545948;


namespace es {


// @cond - Don't document these convenience inline functions or the kernel

// Inverse magnitude of a Vector3d element
static inline double InvMagnitude(Vector3d vec){
    return 1 / vec.norm();
}


// Inverse magnitude of a Vector3d element handling /0 errors to return 0 always
static inline double InvMagnitudeZero(Vector3d vec)
{
    double a = vec.norm();
    if (a == 0.) {
        return 0.;
    } else {
        return 1 / a;
    }
}


// Kernel that applies induction from one vortex line to one control point
static inline Vector3d BiotElementalKernel(Vector3d control_locations, Vector3d start_nodes, Vector3d end_nodes, double gamma, double effective_core_radius_squared)
{

    Vector3d l12 = end_nodes - start_nodes;

    Vector3d r1 = control_locations - start_nodes;

    Vector3d r2 = control_locations - end_nodes;

    Vector3d cross1 = l12.cross(r1);

    // Used invmagnitudezero to solve the problem where the cross product is {0,0,0} due to point p being on the line defined by l12 (r1 collinear with l12)
    Vector3d induction;
    induction = InvMagnitudeZero(cross1) * cross1;
    double mag_l12 = InvMagnitude(l12);
    double mag_r1 = InvMagnitudeZero(r1);
    double mag_r2 = InvMagnitudeZero(r2);
    double costheta_1 = l12.dot(r1) * (mag_l12 * mag_r1);
    double costheta_2 = l12.dot(r2) * (mag_l12 * mag_r2);
    double h = (cross1.norm()) * (mag_l12);
    double magu = gamma * h * (costheta_1 - costheta_2) * invfourpi;
    magu /= sqrt(pow(effective_core_radius_squared, 2) + pow(h, 4));
    induction = induction * magu;

    return induction;

};


// @endcond


Matrix3Xd NaiveBiotSavart(Matrix3Xd startNodes, Matrix3Xd endNodes, Matrix3Xd locations, VectorXd gamma, VectorXd effective_core_radius_squared){

    assert (startNodes.cols() == gamma.size());
    assert (startNodes.cols() == endNodes.cols());
    assert (startNodes.cols() == effective_core_radius_squared.size());
    Index start = 0;
    Index n_control_point = locations.cols();
    Index n_vortex = gamma.size();
    Index n = n_control_point * n_vortex;

    Matrix3Xd induction(3, n_control_point);
    induction.setZero();

    // Parallel version
    tbb::affinity_partitioner partitioner1;
    tbb::affinity_partitioner partitioner2;

    tbb::parallel_for(start, n_control_point, [&](Index control_i){

        Vector3d p = locations.col(control_i);

        Vector3d induction_accumulator;
        induction_accumulator.setZero();
        tbb::spin_mutex induction_mutex;

        tbb::parallel_for(start, n_vortex, [&](Index vortex_i){
            Vector3d A1 = startNodes.col(vortex_i);
            Vector3d B1 = endNodes.col(vortex_i);
            double a = gamma(vortex_i);
            double rc4 = effective_core_radius_squared(vortex_i);
            Vector3d induction_elemental = BiotElementalKernel(p, A1, B1, a, rc4);
            {
                tbb::spin_mutex::scoped_lock lock(induction_mutex);
                induction_accumulator += induction_elemental;
            }
        }, partitioner1);
        induction.block(0, control_i, 3, 1) = induction_accumulator;

    }, partitioner2);

    /* Serial version
    for (long control_i = 0; control_i < n_control_point; control_i++)
    {
        Vector3d p = locations.block(0,control_i,3,1);
        long vortex_i;
        for (vortex_i = 0; vortex_i < n_vortex; vortex_i++) {
            Vector3d A1 = startNodes.block(0, vortex_i, 3, 1);
            Vector3d B1 = endNodes.block(0, vortex_i, 3, 1);
            double a = gamma(vortex_i);
            double rc4 = effective_core_radius_squared(vortex_i);
            induction.block(0, control_i, 3, 1) = induction.block(0, control_i, 3, 1) + BiotElementalKernel(p, A1, B1, a, rc4);
        };

    }
    */

    return induction;
};

};

#endif