corona_equilibrium.h

Go to the documentation of this file.

#pragma once
 
// WARPXM includes
#include "apps/wmapplication.h"
 
// STL includes
#include <cmath>
#include <memory>
 
namespace wxm
{
namespace apps
{
namespace functions
{
namespace mhd
{
class CoronaEquilibrium : public WmApplication
{
public:
    CoronaEquilibrium();
 
    ~CoronaEquilibrium() override;
 
    void setup(const WxCryptSet& wxc) override;
 
    const std::vector<int>& getOutputVariableIndexes(int flag) const override
    {
        return _mhd_fluid;
    }
 
    void evaluate_function(const real* q,
                           const real* aux,
                           const solverVariables_t* pSV,
                           real* result) const override;
 
    void bc_q(const real* q_in,
              const real* aux_in,
              const real* aux_out,
              const solverVariables_t* pFV,
              real* q_out) const override;
 
protected:
    std::vector<int> _mhd_fluid;
 
    std::string _coordinate_system;
    std::string _core_type;
    int _ir, _ith, _iz; // r, theta, and z indices
 
    real _gas_gamma; // gamma
    real _Ai;        // mass
    real _Zi;        // charge
    real _T;         // uniform temperature
    real _a;         // pinch radius
    real _Ib;        // current at radius b
    real _b;         // radius b
    real _wct;       // (omega_c tau)
 
    real _pert;   // perturbation fraction on pressure
    real _lambda; // wavelength of perturbation in z-coordinate
 
    // pointer to subapplications
    std::unique_ptr<WmApplication> _vz_app;
    std::unique_ptr<WmApplication> _vr_app;
    std::unique_ptr<WmApplication> _vth_app;
 
    // Viete function for alpha = 5/3
    real viete(real r, real rp) const
    {
        // Handle the case when r=0
        if (r == 0.0)
        {
            return 0.0; // Will be handled separately in the calling function
        }
 
        // Calculate using the algorithm but avoid standard complex library
        // We'll implement the complex arithmetic directly
 
        // Calculate the squared ratio (rp / r)^2
        real ratio = rp / r;
        real ratio_squared = ratio * ratio;
 
        // Compute the argument for arccos: 6 * sqrt(3) / 5 * (r / rp)^2
        real arccos_arg_real = 6.0 * std::sqrt(3.0) / 5.0 * (r / rp) * (r / rp);
 
        // For complex arccos, we need to handle when arccos_arg_real > 1
        // If arccos_arg_real <= 1, we can use the standard arccos
        // If arccos_arg_real > 1, we use the formula:
        // acos(z) = 0 - i*log(z + i*sqrt(z^2 - 1))
 
        real angle_real;
        real angle_imag;
 
        if (arccos_arg_real <= 1.0)
        {
            // Real case
            angle_real = std::acos(arccos_arg_real) / 3.0;
            angle_imag = 0.0;
        }
        else
        {
            // Complex case
            // For z > 1, acos(z) = 0 - i*log(z + sqrt(z^2 - 1))
            real temp = std::sqrt(arccos_arg_real * arccos_arg_real - 1.0);
            angle_real = 0.0;
            angle_imag = -std::log(arccos_arg_real + temp) / 3.0;
        }
 
        // Compute the cosine of the complex angle
        // cos(a + bi) = cos(a)*cosh(b) - i*sin(a)*sinh(b)
        real cos_term_real = std::cos(angle_real) * std::cosh(angle_imag);
        real cos_term_imag = -std::sin(angle_real) * std::sinh(angle_imag);
 
        // Calculate the final result: (rp/r)^2 * (2/sqrt(3)) * cos_term
        real scale = ratio_squared * (2.0 / std::sqrt(3.0));
        real result_real = scale * cos_term_real;
 
        // Return only the real part
        return result_real;
    }
 
private:
    CoronaEquilibrium& operator=(const CoronaEquilibrium& var);
    CoronaEquilibrium(const CoronaEquilibrium& var);
};
 
} // namespace mhd
} // namespace functions
} // namespace apps
} // namespace wxm