make_field.h

Go to the documentation of this file.

#include "magnetic_field.h"
#include "solovev.h"
#include "guenter.h"
#include "polynomial.h"
#include "toroidal.h"
#include "sheath.h"
#include "modified.h"
#include <dg/file/json_utilities.h>
 
namespace dg{
namespace geo{
 
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);
        }
    }
}
 
 
//TODO The place where these functions are somewhat tested is geometry_diag.cpp Formalise?
 
namespace detail{
 
inline 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));
}
 
inline void createModifiedField(
        const dg::geo::TokamakMagneticField& mag,
        dg::file::WrappedJsonValue jsmod,
        modifier& mod,
        CylindricalFunctorsLvl2& mod_psip,
        CylindricalFunctor& wall, CylindricalFunctor& transition)
{
    std::string m = jsmod.get( "type", "heaviside" ).asString();
    mod = modifier::heaviside;
    description desc = mag.params().getDescription();
    try{
        mod = str2modifier.at( m);
    }catch ( std::out_of_range& err)
    {
        std::string message = "ERROR: Key \"" + m
            + "\" not valid in field:\n\t"
            + jsmod.access_string() + "\"type\" \n";
        throw std::out_of_range(message);
    }
    double RX1 = 0., ZX1 = 0., RX2 = 0., ZX2 = 0.;
    double RO=mag.R0(), ZO=0.;
    if( desc != description::none and desc != description::centeredX)
        dg::geo::findOpoint( mag.get_psip(), RO, ZO);
    if ( desc == description::standardX or desc == description::doubleX)
    {
        // Find first X-point
        RX1 = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
        ZX1 = -1.1*mag.params().elongation()*mag.params().a();
        dg::geo::findXpoint( mag.get_psip(), RX1, ZX1);
    }
    if ( desc == description::doubleX)
    {
        // Find second X-point
        RX2 = mag.R0()-1.1*mag.params().triangularity()*mag.params().a();
        ZX2 = +1.1*mag.params().elongation()*mag.params().a();
        dg::geo::findXpoint( mag.get_psip(), RX2, ZX2);
    }
    switch (mod) {
        default: //none
        {
            wall = mod::DampingRegion( mod::nowhere, mag.psip(), 0, 1, -1);
            transition = mod::MagneticTransition( mod::nowhere, mag.psip(),
                    0,  1, -1);
            mod_psip = mag.get_psip();
            break;
        }
        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(); //given in rho_p
            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){
                // TODO Maybe this can be unified with modifier::sol_pfr_2X
                case description::standardX:
                {
                    double sign0 = +1., sign1 = -1.;
                    detail::transform_psi( mag, psi0, alpha0, sign0);
                    detail::transform_psi( mag, psi1, alpha1, sign1);
                    CylindricalFunctorsLvl2 mod0_psip;
                    mod0_psip = mod::createPsip(
                            mod::everywhere, mag.get_psip(), psi0, alpha0, sign0);
                    mod_psip = mod::createPsip(
                            mod::HeavisideZ( ZX1, -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(ZX1, -1), mag.psip(), psi1, alpha1, -sign1);
                    CylindricalFunctor transition1 = mod::MagneticTransition(
                        mod::HeavisideZ(ZX1, -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);
                    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;
                }
            }
            break;
        }
        case modifier::sol_pfr_2X:
        {
            if( desc != description::doubleX)
                throw Error( Message(_ping_) << "Description must be doubleX");
            unsigned num = 4;
            std::vector<double> psi(num), alpha(num);
            for( unsigned u=0; u<num; u++)
            {
                psi[u] = jsmod["boundary"].get( u, 1.0) .asDouble();
                alpha[u] = jsmod["alpha"].get( u, 0.1) .asDouble();
            }
            std::vector<double> sign = {+1,-1.,+1.,-1};
            for( unsigned u=0; u<num; u++)
                detail::transform_psi( mag, psi[u], alpha[u], sign[u]);
 
            std::vector<std::function< bool( double,double)>> mods = {
                mod::RightSideOf( {RX1,ZX1},{RO,ZO},{RX2,ZX2}),
                mod::Above( {RX1, ZX1}, {RO,ZO}, false), // = Below
                mod::RightSideOf( {RX2,ZX2},{RO,ZO},{RX1,ZX1}),
                mod::Above( {RX2, ZX2}, {RO, ZO}, false) };
 
            CylindricalFunctorsLvl2 mod_psips[num];
            mod_psips[0] = mod::createPsip( mods[0], mag.get_psip(), psi[0],
                alpha[0], sign[0]);
            for( unsigned u=1; u<num; u++)
                mod_psips[u] = mod::createPsip( mods[u], mod_psips[u-1], psi[u],
                    alpha[u], sign[u]);
            mod_psip = mod_psips[num-1];
 
            CylindricalFunctor walls[num], transitions[num];
            for( unsigned u=0; u<4; u++)
            {
                walls[u] = mod::DampingRegion( mods[u], mag.psip(), psi[u],
                    alpha[u], -sign[u]);
                transitions[u] = mod::MagneticTransition( mods[u], mag.psip(),
                    psi[u], alpha[u], sign[u]);
            }
            transition = mod::SetUnion( transitions[0], transitions[1]);
            wall       = mod::SetUnion( walls[0], walls[1]);
            for( unsigned u=2; u<4; u++)
            {
                transition = mod::SetUnion( transition, transitions[u]);
                wall       = mod::SetUnion( wall, walls[u]);
            }
            break;
        }
    }
}
}//namespace detail
 
 
inline TokamakMagneticField createModifiedField(
        dg::file::WrappedJsonValue gs, dg::file::WrappedJsonValue jsmod,
        CylindricalFunctor& wall, CylindricalFunctor& transition)
{
    dg::geo::CylindricalFunctorsLvl2 mod_psip;
    TokamakMagneticField mag = createMagneticField(gs);
    const MagneticFieldParameters& inp = mag.params();
 
    modifier mod;
    detail::createModifiedField( mag, jsmod, mod, mod_psip, wall, transition);
    description desc = inp.getDescription();
    equilibrium equi = inp.getEquilibrium();
    MagneticFieldParameters mod_params{ inp.a(), inp.elongation(),
        inp.triangularity(), equi, mod, desc};
    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);
        }
    }
}
 
 
inline CylindricalFunctor createWallRegion( dg::geo::TokamakMagneticField mag,
    dg::file::WrappedJsonValue jsmod)
{
    CylindricalFunctor wall, transition;
    CylindricalFunctorsLvl2 mod_psip;
    modifier mod;
    detail::createModifiedField( mag, jsmod, mod, mod_psip, wall, transition);
    return wall;
}
inline CylindricalFunctor createWallRegion( dg::file::WrappedJsonValue gs,
    dg::file::WrappedJsonValue jsmod)
{
    return createWallRegion( createMagneticField(gs), jsmod);
}
 
 
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", mod::SOLRegion(mag,wall));
    CylindricalFunctor distP = dg::geo::WallFieldlineDistance( dg::geo::createBHat(
            mag), sheath_walls, (+boundary+1e-3)*2.0*M_PI,
            1e-6, "phi", mod::SOLRegion(mag,wall));
    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