2d2v_kh_instability.h

Go to the documentation of this file.

#ifndef kinetics2d2v_distribution_kelvin_helmholtz_instability_h
#define kinetics2d2v_distribution_kelvin_helmholtz_instability_h
 
// WARPXM includes
#include "apps/wmapplication.h"
 
// STL includes
#include <cmath>
 
namespace wxm
{
namespace apps
{
namespace functions
{
namespace kinetics
{
 
class Kinetics2d2vKelvinHelmholtzInstability : public WmApplication
{
public:
    Kinetics2d2vKelvinHelmholtzInstability();
 
    ~Kinetics2d2vKelvinHelmholtzInstability();
 
    void setup(const WxCryptSet& wxc) override;
 
    const std::vector<int>& getOutputVariableIndexes(int flag) const override
    {
        return _distribution;
    }
 
    void evaluate_function(const real* q,
                           const real* aux,
                           const solverVariables_t* pSV,
                           real* result) const override;
 
    void bc_q_kinetic(const real* q_in,
                      const real* aux_in,
                      const solverVariables_t* pFV,
                      real* q_out) const override;
 
    static real compute_spline_fit(const real x,
                                   const std::vector<real>& xvec,
                                   const std::vector<real>& a,
                                   const std::vector<real>& b,
                                   std::map<real, real>& cache)
    {
        if (cache.count(x)) {
            return cache[x];
        }
 
        size_t n = xvec.size();
        if ((a.size() != n) || (b.size() != 2))
        {
            throw WxExcept(
                "wxm::apps::functions::kinetics::Kinetics2d2vKelvinHelmholtzInstability::"
                "compute_spline_fit() : Vector Size Error!\n");
        }
 
        real val = 0.0;
        for (int i = 0; i < n; i++)
        {
            val += a[i] * std::pow(std::fabs(x - xvec[i]), 3);
        }
        val += b[0] + b[1] * x;
 
        cache[x] = val;
 
        return val;
    }
 
protected:
    // Helper Functions
    real compute_phi_star(const real x) const
    {
        return -_Ex0 * _d_width / 2.0 * std::log(1.0 + std::exp(2.0 * x / _d_width));
    }
 
    real compute_phi_star_double_prime(const real x) const
    {
        return -_Ex0 / _d_width / (1.0 + std::cosh(2.0 * x / _d_width));
    }
 
    real compute_gi(const real x) const
    {
        return std::pow(1.0 + std::exp(2.0 * x / _d_width), _b_gradient * _d_width / 2.0);
    }
 
    real compute_ge(const real x) const
    {
        const real ni_aux_fit = compute_spline_fit(
            x, _ni_aux_fit_xvec, _ni_aux_fit_acoeffs, _ni_aux_fit_bcoeffs, _ni_aux_fit_cache);
 
        return -_Zi / _Ze * ni_aux_fit - 1.0 / _Ze * _omega_c_tau / _omega_p_tau /
                                             _omega_p_tau *
                                             compute_phi_star_double_prime(x);
    }
 
    void compute_fi(const solverVariables_t* pSV, real* result) const;
 
    void compute_fe(const solverVariables_t* pSV, real* result) const;
 
    // Variables
    std::vector<int> _distribution; // output variable
 
    std::string _species;      // "ions" or "electrons"
    real _Ex0;                 // E0 / omega_c_tau
    real _d_width;             // shear layer width
    real _b_gradient;          // density gradient term
    real _temperature_uniform; // T
    real _Zi;                  // ion charge
    real _Ze;                  // electron charge, should only ever be -1
    real _Ai;                  // ion mass
    real _Ae;                  // electron mass
    real _Omega_c_i;           // omega_c_tau * Zi / Ai
    real _Omega_c_e;           // omega_c_tau * Ze / Ae
    real _omega_p_tau;         // (omega_p tau)
    real _omega_c_tau;         // (omega_c tau)
 
    std::vector<real> _perturbation_center;
    real _perturbation_wavelength;
    real _perturbation_amplitude;
 
    std::vector<real> _ni_aux_fit_xvec;
    std::vector<real> _ni_aux_fit_acoeffs;
    std::vector<real> _ni_aux_fit_bcoeffs;
 
    std::vector<real> _phi_fit_xvec;
    std::vector<real> _phi_fit_acoeffs;
    std::vector<real> _phi_fit_bcoeffs;
 
    std::vector<real> _Ay_fit_xvec;
    std::vector<real> _Ay_fit_acoeffs;
    std::vector<real> _Ay_fit_bcoeffs;
 
    mutable std::map<real, real> _Ay_fit_cache;
    mutable std::map<real, real> _ni_aux_fit_cache;
    mutable std::map<real, real> _phi_fit_cache;
 
private:
    Kinetics2d2vKelvinHelmholtzInstability&
    operator=(const Kinetics2d2vKelvinHelmholtzInstability& var);
    Kinetics2d2vKelvinHelmholtzInstability(
        const Kinetics2d2vKelvinHelmholtzInstability& var);
};
 
} // namespace kinetics
} // namespace functions
} // namespace apps
} // namespace wxm
 
#endif