make_field.h

Go to the documentation of this file.

#ifdef JSONCPP_VERSION_STRING
#include "magnetic_field.h"
#include "solovev.h"
#include "guenter.h"
#include "polynomial.h"
#include "toroidal.h"
#include "fieldaligned.h"
#include <dg/file/json_utilities.h>
 
namespace dg{
namespace geo{
 
static inline TokamakMagneticField createMagneticField( dg::file::WrappedJsonValue gs)
{
    std::string e = gs.get( "equilibrium", "solovev" ).asString();
    equilibrium equi = equilibrium::solovev;
    try{
        equi = str2equilibrium.at( e);
    }catch ( std::out_of_range& err)
    {
        std::string message = "ERROR: Key \"" + e
            + "\" not valid in field:\n\t"
            + gs.access_string() + "\"equilibrium\" \n";
        throw std::out_of_range(message);
    }
    switch( equi){
        case equilibrium::polynomial:
        {
            polynomial::Parameters gp( gs);
            return createPolynomialField( gp);
        }
        case equilibrium::toroidal:
        {
            double R0 = gs.get( "R_0", 10.0).asDouble();
            return createToroidalField( R0);
        }
        case equilibrium::guenter:
        {
            double I0 = gs.get( "I_0", 20.0).asDouble();
            double R0 = gs.get( "R_0", 10.0).asDouble();
            return createGuenterField( R0, I0);
        }
        case equilibrium::circular:
        {
            double I0 = gs.get( "I_0", 20.0).asDouble();
            double R0 = gs.get( "R_0", 10.0).asDouble();
            double a = gs.get( "a", 1.0).asDouble();
            double b = gs.get( "b", 1.0).asDouble();
            return createCircularField( R0, I0, a, b);
        }
#ifdef BOOST_VERSION
        case equilibrium::taylor:
        {
            solovev::Parameters gp( gs);
            return createTaylorField( gp);
        }
#endif
        default:
        {
            solovev::Parameters gp( gs);
            return createSolovevField( gp);
        }
    }
}
 
 
namespace detail{
void transform_psi( TokamakMagneticField mag, double& psi0, double& alpha0, double& sign0)
{
    double RO=mag.R0(), ZO=0.;
    dg::geo::findOpoint( mag.get_psip(), RO, ZO);
    double psipO = mag.psip()( RO, ZO);
    double wall_psi0p = (1.-psi0*psi0)*psipO;
    double wall_alpha0p = -(2.*psi0+alpha0)*alpha0*psipO;
    psi0 = wall_psi0p + sign0*wall_alpha0p/2.;
    alpha0 = fabs( wall_alpha0p/2.);
    sign0 = sign0*((psipO>0)-(psipO<0));
}
}//namespace detail
 
 
static inline TokamakMagneticField createModifiedField(
        dg::file::WrappedJsonValue gs, dg::file::WrappedJsonValue jsmod,
        CylindricalFunctor& wall, CylindricalFunctor& transition)
{
    TokamakMagneticField mag = createMagneticField( gs);
    const MagneticFieldParameters& inp = mag.params();
    description desc = inp.getDescription();
    equilibrium equi = inp.getEquilibrium();
    std::string m = jsmod.get( "type", "heaviside" ).asString();
    modifier mod = str2modifier.at( m);
    MagneticFieldParameters mod_params{ inp.a(), inp.elongation(),
        inp.triangularity(), equi, mod, desc};
    CylindricalFunctorsLvl2 mod_psip;
    switch (mod) {
        default: //none
        {
            wall = mod::DampingRegion( mod::nowhere, mag.psip(), 0, 1, -1);
            transition = mod::MagneticTransition( mod::nowhere, mag.psip(),
                    0,  1, -1);
            return mag;
        }
        case modifier::heaviside:
        {
            double psi0 = jsmod.get( "boundary", 1.1 ).asDouble();
            double alpha = jsmod.get( "alpha", 0.2 ).asDouble();
            double sign = +1;
            if( desc == description::standardX || desc == description::standardO ||
                    desc == description::doubleX)
                detail::transform_psi( mag, psi0, alpha, sign);
            else
                sign = jsmod.get( "sign", -1. ).asDouble();
 
            mod_psip = mod::createPsip( mod::everywhere, mag.get_psip(), psi0, alpha, sign);
            wall = mod::DampingRegion( mod::everywhere, mag.psip(), psi0, alpha, -sign);
            transition = mod::MagneticTransition( mod::everywhere, mag.psip(), psi0, alpha, sign);
            break;
        }
        case modifier::sol_pfr:
        {
            double psi0     = jsmod["boundary"].get( 0, 1.1 ).asDouble();
            double alpha0   = jsmod["alpha"].get(0, 0.2 ).asDouble();
            double psi1     = jsmod["boundary"].get(1, 0.97 ).asDouble();
            double alpha1   = jsmod[ "alpha"].get(1, 0.2 ).asDouble();
            switch( desc){
                case description::standardX:
                {
                    double sign0 = +1., sign1 = -1.;
                    detail::transform_psi( mag, psi0, alpha0, sign0);
                    detail::transform_psi( mag, psi1, alpha1, sign1);
                    //we can find the X-point
                    double RX = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
                    double ZX = -1.1*mag.params().elongation()*mag.params().a();
                    dg::geo::findXpoint( mag.get_psip(), RX, ZX);
                    CylindricalFunctorsLvl2 mod0_psip;
                    mod0_psip = mod::createPsip(
                            mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
                    mod_psip = mod::createPsip(
                            mod::HeavisideZ( ZX, -1), mod0_psip, psi1, alpha1, sign1);
                    CylindricalFunctor wall0 = mod::DampingRegion( mod::everywhere, mag.psip(), psi0, alpha0, -sign0);
                    CylindricalFunctor transition0 = mod::MagneticTransition( mod::everywhere, mag.psip(), psi0, alpha0, sign0);
                    CylindricalFunctor wall1 = mod::DampingRegion( mod::HeavisideZ(ZX, -1), mag.psip(), psi1, alpha1, -sign1);
                    CylindricalFunctor transition1 = mod::MagneticTransition( mod::HeavisideZ(ZX, -1), mag.psip(), psi1, alpha1, sign1);
                    wall = mod::SetUnion( wall0, wall1);
                    transition = mod::SetUnion( transition0, transition1);
                    break;
                }
                case description::doubleX:
                {
                    double sign0 = +1., sign1 = -1.;
                    detail::transform_psi( mag, psi0, alpha0, sign0);
                    detail::transform_psi( mag, psi1, alpha1, sign1);
                    //we can find the X-point
                    double RX1 = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
                    double ZX1 = -1.1*mag.params().elongation()*mag.params().a();
                    dg::geo::findXpoint( mag.get_psip(), RX1, ZX1);
                    double RX2 = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
                    double ZX2 = +1.1*mag.params().elongation()*mag.params().a();
                    dg::geo::findXpoint( mag.get_psip(), RX2, ZX2);
                    CylindricalFunctorsLvl2 mod0_psip, mod1_psip;
                    mod0_psip = mod::createPsip(
                            mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
                    mod1_psip = mod::createPsip(
                            mod::HeavisideZ( ZX1, -1), mod0_psip, psi1, alpha1, sign1);
                    mod_psip = mod::createPsip(
                            mod::HeavisideZ( ZX2, +1), mod1_psip, psi1, alpha1, sign1);
                    CylindricalFunctor wall0 = mod::DampingRegion( mod::everywhere, mag.psip(), psi0, alpha0, -sign0);
                    CylindricalFunctor wall1 = mod::DampingRegion( mod::HeavisideZ(ZX1, -1), mag.psip(), psi1, alpha1, -sign1);
                    CylindricalFunctor wall2 = mod::DampingRegion( mod::HeavisideZ(ZX2, +1), mag.psip(), psi1, alpha1, -sign1);
                    CylindricalFunctor transition0 = mod::MagneticTransition( mod::everywhere, mag.psip(), psi0, alpha0, sign0);
                    CylindricalFunctor transition1 = mod::MagneticTransition( mod::HeavisideZ(ZX1, -1), mag.psip(), psi1, alpha1, sign1);
                    CylindricalFunctor transition2 = mod::MagneticTransition( mod::HeavisideZ(ZX2, +1), mag.psip(), psi1, alpha1, sign1);
                    transition = mod::SetUnion( mod::SetUnion( transition0, transition1), transition2);
                    wall = mod::SetUnion( mod::SetUnion( wall0, wall1), wall2);
                    break;
                }
                default:
                {
                    double sign0 = jsmod[ "sign"].get( 0, -1. ).asDouble();
                    double sign1 = jsmod[ "sign"].get( 1, +1. ).asDouble();
                    CylindricalFunctorsLvl2 mod0_psip;
                    mod0_psip = mod::createPsip(
                            mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
                    mod_psip = mod::createPsip(
                            mod::everywhere, mod0_psip, psi1, alpha1, sign1);
                    CylindricalFunctor wall0 = mod::DampingRegion( mod::everywhere, mag.psip(), psi0, alpha0, sign0);
                    CylindricalFunctor transition0 = mod::MagneticTransition( mod::everywhere, mag.psip(), psi0, alpha0, sign0);
                    CylindricalFunctor wall1 = mod::DampingRegion( mod::everywhere, mag.psip(), psi1, alpha1, sign1);
                    CylindricalFunctor transition1 = mod::MagneticTransition( mod::everywhere, mag.psip(), psi1, alpha1, sign1);
                    wall = mod::SetUnion( wall0, wall1);
                    transition = mod::SetUnion( transition0, transition1);
                    break;
                }
            }
        }
    }
    switch( equi){
        case equilibrium::solovev:
        {
            solovev::Parameters gp( gs);
            return TokamakMagneticField( gp.R_0,mod_psip,
                solovev::createIpol( gp, mod_psip), mod_params);
        }
        default:
        {
            return TokamakMagneticField( mag.R0(), mod_psip,
                    mag.get_ipol(), mod_params);
        }
    }
}
 
 
static inline CylindricalFunctor createWallRegion( dg::file::WrappedJsonValue gs, dg::file::WrappedJsonValue jsmod)
{
    CylindricalFunctor wall, transition;
    TokamakMagneticField mag = createModifiedField( gs, jsmod, wall, transition);
    return wall;
}
 
static inline void createSheathRegion(
    dg::file::WrappedJsonValue jsmod, TokamakMagneticField mag,
    CylindricalFunctor wall, dg::Grid2d& sheath_walls,
    CylindricalFunctor& sheath)
{
    double R0 = sheath_walls.x0(), R1 = sheath_walls.x1();
    double Z0 = sheath_walls.y0(), Z1 = sheath_walls.y1();
    Grid1d gR1d ( R0, R1, 1, 100);
    Grid1d gZ1d ( Z0, Z1, 1, 100);
    std::array<bool,2> sheathR = {false,false}, sheathZ = {false,false};
    for ( unsigned i=0; i<100; i++)
    {
        if( wall( R0+i*gR1d.h(), Z0) == 0)
            sheathR[0] = true;
        if( wall( R0+i*gR1d.h(), Z1) == 0)
            sheathR[1] = true;
        if( wall( R0, Z0 + i*gZ1d.h()) == 0)
            sheathZ[0] = true;
        if( wall( R1, Z0 + i*gZ1d.h()) == 0)
            sheathZ[1] = true;
    }
    if( false == sheathR[0]) Z0 = -1e10;
    if( false == sheathR[1]) Z1 = 1e10;
    if( false == sheathZ[0]) R0 = -1e10;
    if( false == sheathZ[1]) R1 = 1e10;
    sheath_walls = dg::Grid2d( R0, R1, Z0, Z1, 1,1,1);
    double boundary = jsmod.get( "boundary", 0.0625 ).asDouble(); // 1/16
    double alpha    = jsmod.get( "alpha", 0.015625 ).asDouble(); // 1/64
    CylindricalFunctor distM = dg::geo::WallFieldlineDistance( dg::geo::createBHat(
            mag), sheath_walls, (-boundary-1e-3)*2.0*M_PI,
            1e-6, "phi");
    CylindricalFunctor distP = dg::geo::WallFieldlineDistance( dg::geo::createBHat(
            mag), sheath_walls, (+boundary+1e-3)*2.0*M_PI,
            1e-6, "phi");
    dg::PolynomialHeaviside polyM( -boundary*2.*M_PI + alpha*M_PI, alpha*M_PI, +1);
    dg::PolynomialHeaviside polyP(  boundary*2.*M_PI - alpha*M_PI, alpha*M_PI, -1);
    auto sheathM = dg::compose( polyM, distM); //positive (because distance)
    auto sheathP = dg::compose( polyP, distP);
    sheath = mod::SetUnion( sheathM, sheathP);
    sheath = mod::SetIntersection( mod::SetNot( wall), sheath);
}
 
} //namespace geo
}//namespace dg
#endif //JSONCPP_VERSION_STRING