kernel_functions.h

// This file is part of openPSTD.                                       //
//                                                                      //
// openPSTD is free software: you can redistribute it and/or modify     //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or    //
// (at your option) any later version.                                  //
//                                                                      //
// openPSTD is distributed in the hope that it will be useful,          //
// but WITHOUT ANY WARRANTY; without even the implied warranty of       //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        //
// GNU General Public License for more details.                         //
//                                                                      //
// You should have received a copy of the GNU General Public License    //
// along with openPSTD.  If not, see <http://www.gnu.org/licenses/>.    //
//                                                                      //

//
// Date:
//      3-9-2015
//
// Authors:
//      Omar Richardson
//      Louis van Harten
//
// Purpose:
//      Contains procedural functions that don't fit with any one class.
//


#ifndef OPENPSTD_KERNEL_FUNCTIONS_H
#define OPENPSTD_KERNEL_FUNCTIONS_H

#include <Eigen/Dense>
#include <complex>
#include <fftw3.h>
#include <iostream>
#include <vector>
#include <map>
#include <math.h>
#include <algorithm>
#include "../KernelInterface.h"
#include "Geometry.h"

namespace OpenPSTD {
    namespace Kernel {

        template<typename T>
        std::vector<T> arange(T start, T stop, T step = 1) {
            std::vector<T> list;
            for (T i = start; i < stop; i += step)
                list.push_back(i);
            return list;
        }

        const float EPSILON = 1E-5; // Depends on float/double


        enum class CalcDirection {
            X = 0, Y = 1
        };

        enum class CalculationType {
            PRESSURE, VELOCITY
        };

        enum class Direction {
            LEFT, RIGHT, TOP, BOTTOM
        };

        const std::vector<Direction> all_directions = {Direction::LEFT, Direction::RIGHT, Direction::TOP,
                                                       Direction::BOTTOM};

        const std::vector<CalcDirection> all_calc_directions = {CalcDirection::X, CalcDirection::Y};
        const std::vector<CalculationType> all_calculation_types = {CalculationType::PRESSURE,
                                                                    CalculationType::VELOCITY};

        const std::vector<float> rk_coefficients = {8.91421261e-4f, 7555704391e-3f, 4.0919732041e-2f,
                                                    1.65919771368e-1f, 5e-1, 1.f}; // Temporary until bugfix


        Direction get_opposite(Direction direction);

        CalcDirection get_orthogonal(CalcDirection direction);

        CalcDirection direction_to_calc_direction(Direction direction);

        struct RhoArray {
            Eigen::ArrayXXf pressure;
            Eigen::ArrayXXf velocity;
        };

        Eigen::ArrayXXf spatderp3(Eigen::ArrayXXf p1, Eigen::ArrayXXf p2,
                                  Eigen::ArrayXXf p3, Eigen::ArrayXcf derfact,
                                  RhoArray rho_array, Eigen::ArrayXf window, int wlen,
                                  CalculationType ct, CalcDirection direct);

        Eigen::ArrayXXf spatderp3(Eigen::ArrayXXf p1, Eigen::ArrayXXf p2,
                                  Eigen::ArrayXXf p3, Eigen::ArrayXcf derfact,
                                  RhoArray rho_array, Eigen::ArrayXf window, int wlen,
                                  CalculationType ct, CalcDirection direct,
                                  fftwf_plan plan, fftwf_plan plan_inv);

        RhoArray get_rho_array(const float rho1, const float rho_self, const float rho2);

        float get_grid_spacing(PSTDSettings cnf);

        Eigen::ArrayXf get_window_coefficients(int window_size, int patch_error);

        int next_2_power(float n);

        bool is_approx(float a, float b);

        void debug(std::string msg);

    }
}
#endif //OPENPSTD_KERNEL_FUNCTIONS_H