solovev.h

Go to the documentation of this file.

#pragma once
 
#include <iostream>
#include <cmath>
#include <vector>
 
#include "dg/blas.h"
 
#include "dg/algorithm.h"
#include "solovev_parameters.h"
#include "magnetic_field.h"
#include "modified.h"
 
 
namespace dg
{
namespace geo
{
namespace solovev
{
 
struct Psip: public aCylindricalFunctor<Psip>
{
    Psip( const Parameters& gp ): m_R0(gp.R_0), mA(gp.A), m_pp(gp.pp), mc(gp.c) {
        m_prev = std::make_shared< std::array<double,3>>(std::array<double,3>{ 0,0,0});
    }
    double do_compute(double R, double Z) const
    {
        // Optimization rationale: The way we compute magnetic field terms
        // through the TokamakMagneticField class and e.g. the BHatR class
        // leads to repeated evaluations of Psip, PsipR, etc. at the same
        // point. We thus let the Psip classes remember the result of the
        // previous call to avoid recomputing the same point. Since the
        // different calls are made through different copies of Psip we need to
        // store the previous results in a shared_ptr such that all copies of
        // Psip have access to it
        if( R == (*m_prev)[0] && Z == (*m_prev)[1])
            return (*m_prev)[2];
 
        double Rn = R / m_R0, Rn2 = Rn * Rn, lgRn = log(Rn);
        double Zn = Z/m_R0;
        // Copied from Mathematica ...
        (*m_prev)[2] =   m_R0*m_pp*(mc[0] + Rn2 * ((lgRn * mA) / 2. + mc[1] +
            Zn * (mc[8] + Zn *
                (-4 * mc[3] - 9 * mc[4] +
                    Zn * (Zn * (8 * mc[5] - 140 * mc[6]) -
                        4 * mc[10]))) +
            lgRn * (-mc[2] +
                Zn * (-3 * mc[9] +
                    Zn * (-12 * mc[4] +
                        Zn * (-120 * Zn * mc[6] - 80 * mc[11]))
                    )) + Rn2 *
            (0.125 - mA / 8. + mc[3] +
                Rn2 * (mc[5] - 15 * lgRn * mc[6]) +
                Zn * (Zn * (-12 * mc[5] + 75 * mc[6]) +
                    3 * mc[10] - 45 * mc[11]) +
                lgRn * (3 * mc[4] +
                    Zn * (180 * Zn * mc[6] + 60 * mc[11]))))\
            + Zn * (mc[7] + Zn *
                (mc[2] + Zn *
                    (mc[9] +
                        Zn * (2 * mc[4] +
                            Zn * (8 * Zn * mc[6] + 8 * mc[11])))))
                      );
        (*m_prev)[0] = R, (*m_prev)[1] = Z;
        return (*m_prev)[2];
    }
  private:
    double m_R0, mA, m_pp;
    std::vector<double> mc;
    std::shared_ptr<std::array<double,3>> m_prev;
};
 
struct PsipR: public aCylindricalFunctor<PsipR>
{
    PsipR( const Parameters& gp ): m_R0(gp.R_0), mA(gp.A), m_pp(gp.pp), mc(gp.c) {
        m_prev = std::make_shared< std::array<double,3>>(std::array<double,3>{ 0,0,0});
    }
    double do_compute(double R, double Z) const
    {
        if( R == (*m_prev)[0] && Z == (*m_prev)[1])
            return (*m_prev)[2];
        double Rn = R / m_R0, Rn2 = Rn * Rn, lgRn = log(Rn);
        double Zn = Z / m_R0;
        // Copied from Mathematica ...
        (*m_prev)[2] = m_pp * (mA * Rn * (0.5 + lgRn - Rn2 / 2.) +
            Rn * (2 * mc[1] - mc[2] -
                8 * Zn * Zn * mc[3] +
                lgRn * (-2 * mc[2] +
                    Zn * (-6 * mc[9] +
                        Zn * (-24 * mc[4] +
                            Zn * (-240 * Zn * mc[6] - 160 * mc[11])
                            ))) +
                Zn * (2 * mc[8] - 3 * mc[9] +
                    Zn * (-30 * mc[4] +
                        Zn * (Zn * (16 * mc[5] - 400 * mc[6]) -
                            8 * mc[10] - 80 * mc[11]))) +
                Rn2 * (0.5 + 4 * mc[3] + 3 * mc[4] +
                    Rn2 * (6 * mc[5] - 15 * mc[6] -
                        90 * lgRn * mc[6]) +
                    Zn * (Zn * (-48 * mc[5] + 480 * mc[6]) +
                        12 * mc[10] - 120 * mc[11]) +
                    lgRn * (12 * mc[4] +
                        Zn * (720 * Zn * mc[6] + 240 * mc[11]))))
            );
        (*m_prev)[0] = R, (*m_prev)[1] = Z;
        return (*m_prev)[2];
    }
  private:
    double m_R0, mA, m_pp;
    std::vector<double> mc;
    std::shared_ptr<std::array<double,3>> m_prev;
};
struct PsipRR: public aCylindricalFunctor<PsipRR>
{
    PsipRR( const Parameters& gp ): m_R0(gp.R_0), mA(gp.A), m_pp(gp.pp), mc(gp.c) {
        m_prev = std::make_shared< std::array<double,3>>(std::array<double,3>{ 0,0,0});
    }
    double do_compute(double R, double Z) const
    {
        if( R == (*m_prev)[0] && Z == (*m_prev)[1])
            return (*m_prev)[2];
        double Rn = R / m_R0, Rn2 = Rn * Rn, lgRn = log(Rn);
        double Zn = Z / m_R0;
        // Copied from Mathematica ...
        (*m_prev)[2] = m_pp/m_R0 * (mA * (1.5 + lgRn - (3 * Rn2) / 2.) + 2 * mc[1] -
            3 * mc[2] - 8 * Zn * Zn * mc[3] +
            lgRn * (-2 * mc[2] +
                Zn * (-6 * mc[9] +
                    Zn * (-24 * mc[4] +
                        Zn * (-240 * Zn * mc[6] - 160 * mc[11]))))
            + Zn * (2 * mc[8] - 9 * mc[9] +
                Zn * (-54 * mc[4] +
                    Zn * (Zn * (16 * mc[5] - 640 * mc[6]) -
                        8 * mc[10] - 240 * mc[11]))) +
            Rn2 * (1.5 + 12 * mc[3] + 21 * mc[4] +
                Rn2 * (30 * mc[5] - 165 * mc[6] -
                    450 * lgRn * mc[6]) +
                Zn * (Zn * (-144 * mc[5] + 2160 * mc[6]) +
                    36 * mc[10] - 120 * mc[11]) +
                lgRn * (36 * mc[4] +
                    Zn * (2160 * Zn * mc[6] + 720 * mc[11])))
            );
        (*m_prev)[0] = R, (*m_prev)[1] = Z;
        return (*m_prev)[2];
    }
  private:
    double m_R0, mA, m_pp;
    std::vector<double> mc;
    std::shared_ptr<std::array<double,3>> m_prev;
};
struct PsipZ: public aCylindricalFunctor<PsipZ>
{
    PsipZ( const Parameters& gp ): m_R0(gp.R_0), m_pp(gp.pp), mc(gp.c) {
        m_prev = std::make_shared< std::array<double,3>>(std::array<double,3>{ 0,0,0});
    }
    double do_compute(double R, double Z) const
    {
        if( R == (*m_prev)[0] && Z == (*m_prev)[1])
            return (*m_prev)[2];
        double Rn = R / m_R0, Rn2 = Rn * Rn, lgRn = log(Rn);
        double Zn = Z / m_R0;
        // Copied from Mathematica ...
        (*m_prev)[2] = m_pp * (mc[7] + Rn2 * (mc[8] - 3 * lgRn * mc[9] +
            Rn2 * (3 * mc[10] - 45 * mc[11] +
                60 * lgRn * mc[11])) +
            Zn * (2 * mc[2] + Rn2 *
                (-8 * mc[3] + (-18 - 24 * lgRn) * mc[4] +
                    Rn2 * (-24 * mc[5] + 150 * mc[6] +
                        360 * lgRn * mc[6])) +
                Zn * (3 * mc[9] +
                    Rn2 * (-12 * mc[10] - 240 * lgRn * mc[11]) +
                    Zn * (8 * mc[4] +
                        Rn2 * (32 * mc[5] - 560 * mc[6] -
                            480 * lgRn * mc[6]) +
                        Zn * (48 * Zn * mc[6] + 40 * mc[11]))))
            );
        (*m_prev)[0] = R, (*m_prev)[1] = Z;
        return (*m_prev)[2];
 
    }
  private:
    double m_R0, m_pp;
    std::vector<double> mc;
    std::shared_ptr<std::array<double,3>> m_prev;
};
struct PsipZZ: public aCylindricalFunctor<PsipZZ>
{
    PsipZZ( const Parameters& gp): m_R0(gp.R_0), m_pp(gp.pp), mc(gp.c) {
        m_prev = std::make_shared< std::array<double,3>>(std::array<double,3>{ 0,0,0});
    }
    double do_compute(double R, double Z) const
    {
        if( R == (*m_prev)[0] && Z == (*m_prev)[1])
            return (*m_prev)[2];
        double Rn = R / m_R0, Rn2 = Rn * Rn, lgRn = log(Rn);
        double Zn = Z / m_R0;
        // Copied from Mathematica ...
        (*m_prev)[2] = m_pp/m_R0 * (2 * mc[2] + 24 * Zn * Zn * mc[4] +
            Zn * (6 * mc[9] + Zn * Zn *
                (240 * Zn * mc[6] + 160 * mc[11])) +
            Rn2 * (-8 * mc[3] + (-18 - 24 * lgRn) * mc[4] +
                Rn2 * (-24 * mc[5] + 150 * mc[6] +
                    360 * lgRn * mc[6]) +
                Zn * (Zn * (96 * mc[5] - 1680 * mc[6] -
                    1440 * lgRn * mc[6]) - 24 * mc[10] -
                    480 * lgRn * mc[11]))
            );
        (*m_prev)[0] = R, (*m_prev)[1] = Z;
        return (*m_prev)[2];
    }
  private:
    double m_R0, m_pp;
    std::vector<double> mc;
    std::shared_ptr<std::array<double,3>> m_prev;
};
struct PsipRZ: public aCylindricalFunctor<PsipRZ>
{
    PsipRZ( const Parameters& gp ): m_R0(gp.R_0), m_pp(gp.pp), mc(gp.c) {
        m_prev = std::make_shared< std::array<double,3>>(std::array<double,3>{ 0,0,0});
    }
    double do_compute(double R, double Z) const
    {
        if( R == (*m_prev)[0] && Z == (*m_prev)[1])
            return (*m_prev)[2];
        double Rn = R / m_R0, Rn2 = Rn * Rn, lgRn = log(Rn);
        double Zn = Z / m_R0;
        // Copied from Mathematica ...
        (*m_prev)[2] = m_pp/m_R0 * (Rn * (2 * mc[8] - 3 * mc[9] - 6 * lgRn * mc[9] +
            Rn2 * (Zn * (-96 * mc[5] + 960 * mc[6] +
                1440 * lgRn * mc[6]) + 12 * mc[10] -
                120 * mc[11] + 240 * lgRn * mc[11]) +
            Zn * (-16 * mc[3] + (-60 - 48 * lgRn) * mc[4] +
                Zn * (Zn * (64 * mc[5] - 1600 * mc[6] -
                    960 * lgRn * mc[6]) - 24 * mc[10] -
                    240 * mc[11] - 480 * lgRn * mc[11])))
            );
        (*m_prev)[0] = R, (*m_prev)[1] = Z;
        return (*m_prev)[2];
    }
  private:
    double m_R0, m_pp;
    std::vector<double> mc;
    std::shared_ptr<std::array<double,3>> m_prev;
};
 
struct Ipol: public aCylindricalFunctor<Ipol>
{
    Ipol( const Parameters &gp, const std::function<double(double,double)>& psip ):  m_R0(gp.R_0), m_A(gp.A), m_pp(gp.pp), m_pi(gp.pi), m_psip(psip) {
        if( gp.pp == 0.)
            m_pp = 1.; //safety measure to avoid divide by zero errors
    }
    double do_compute(double R, double Z) const
    {
        return m_pi*sqrt(-2.*m_A* m_psip(R,Z) /m_R0/m_pp + 1.);
    }
  private:
    double m_R0, m_A, m_pp, m_pi;
    std::function<double(double,double)> m_psip;
};
struct IpolR: public aCylindricalFunctor<IpolR>
{
    IpolR( const Parameters& gp, const std::function<double(double,double)>& psip, std::function<double(double,double)> psipR ):
        m_R0(gp.R_0), m_A(gp.A), m_pp(gp.pp), m_pi(gp.pi), m_psip(psip), m_psipR(psipR) {
        if( gp.pp == 0.)
            m_pp = 1.; //safety measure to avoid divide by zero errors
    }
    double do_compute(double R, double Z) const
    {
        return -m_pi/sqrt(-2.*m_A* m_psip(R,Z) /m_R0/m_pp + 1.)*(m_A*m_psipR(R,Z)/m_R0/m_pp);
    }
  private:
    double m_R0, m_A, m_pp, m_pi;
    std::function<double(double,double)> m_psip, m_psipR;
};
struct IpolZ: public aCylindricalFunctor<IpolZ>
{
    IpolZ(  const Parameters& gp, const std::function<double(double,double)>& psip, std::function<double(double,double)> psipZ ):
        m_R0(gp.R_0), m_A(gp.A), m_pp(gp.pp), m_pi(gp.pi), m_psip(psip), m_psipZ(psipZ) {
        if( gp.pp == 0.)
            m_pp = 1.; //safety measure to avoid divide by zero errors
    }
    double do_compute(double R, double Z) const
    {
        return -m_pi/sqrt(-2.*m_A* m_psip(R,Z) /m_R0/m_pp + 1.)*(m_A*m_psipZ(R,Z)/m_R0/m_pp);
    }
  private:
    double m_R0, m_A, m_pp, m_pi;
    std::function<double(double,double)> m_psip, m_psipZ;
};
 
inline dg::geo::CylindricalFunctorsLvl2 createPsip( const Parameters& gp)
{
    return CylindricalFunctorsLvl2( Psip(gp), PsipR(gp), PsipZ(gp),
        PsipRR(gp), PsipRZ(gp), PsipZZ(gp));
}
inline dg::geo::CylindricalFunctorsLvl1 createIpol( const Parameters& gp, const CylindricalFunctorsLvl1& psip)
{
    return CylindricalFunctorsLvl1(
            solovev::Ipol(gp, psip.f()),
            solovev::IpolR(gp,psip.f(), psip.dfx()),
            solovev::IpolZ(gp,psip.f(), psip.dfy()));
}
 
 
} //namespace solovev
 
inline dg::geo::TokamakMagneticField createSolovevField(
    const dg::geo::solovev::Parameters& gp)
{
    MagneticFieldParameters params = { gp.a, gp.elongation, gp.triangularity,
            equilibrium::solovev, modifier::none, str2description.at( gp.description)};
    auto psip = solovev::createPsip(gp); // make sure prev is shared
    return TokamakMagneticField( gp.R_0, psip,
        solovev::createIpol(gp, psip), params);
}
inline dg::geo::TokamakMagneticField createModifiedSolovevField(
    const dg::geo::solovev::Parameters& gp, double psi0, double alpha, double sign = -1)
{
    MagneticFieldParameters params = { gp.a, gp.elongation, gp.triangularity,
            equilibrium::solovev, modifier::heaviside, str2description.at( gp.description)};
    auto psip = solovev::createPsip(gp); // make sure prev is shared
    return TokamakMagneticField( gp.R_0,
            mod::createPsip( mod::everywhere, psip, psi0, alpha, sign),
        solovev::createIpol( gp, mod::createPsip( mod::everywhere, psip, psi0, alpha, sign)),
        params);
}
 
} //namespace geo
} //namespace dg