analytic_3dmag.hpp

Go to the documentation of this file.

#ifndef ANALYTIC_3D_MAG_HPP
#define ANALYTIC_3D_MAG_HPP

#include "globals.hpp"
#include "equil.hpp"

// Provides setup and get psi for an analytic example of psi
// get_analytic_psi() will be replaced with bicubic interp eventually

namespace{

double get_circular_analytic_psi3d(const Equilibrium& equil, double r, double z, double phi){
    double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
    double zp = (z - equil.axis.z)/equil.bounds.get_z_width();

    double phi_factor = 1.0 + 0.1*sin(phi);

    return (rp*rp + zp*zp)*phi_factor;
}

double get_analytic_dpsidr(const Equilibrium& equil, double r, double z, double phi){
    double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
    double zp = (z - equil.axis.z)/equil.bounds.get_z_width();

    double phi_factor = 1.0 + 0.1*sin(phi);

    return 2*rp/equil.bounds.get_r_width()*phi_factor;
}

double get_analytic_dpsidz(const Equilibrium& equil, double r, double z, double phi){
    double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
    double zp = (z - equil.axis.z)/equil.bounds.get_z_width();

    double phi_factor = 1.0 + 0.1*sin(phi);

    return 2*zp/equil.bounds.get_z_width()*phi_factor;
}

double get_analytic_dpsidphi(const Equilibrium& equil, double r, double z, double phi){
    double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
    double zp = (z - equil.axis.z)/equil.bounds.get_z_width();

    return (rp*rp + zp*zp)*0.1*cos(phi);
}

// Br = -dpsi/dz / r
double get_analytic_Br(const Equilibrium& equil, double r, double z, double phi){
    double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
    double zp = (z - equil.axis.z)/equil.bounds.get_z_width();

    double phi_factor = 1.0 + 0.1*sin(phi);

    return -2*zp/equil.bounds.get_z_width()*phi_factor/r;
}

// Bz = dpsi/dr / r
double get_analytic_Bz(const Equilibrium& equil, double r, double z, double phi){
    double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
    double zp = (z - equil.axis.z)/equil.bounds.get_z_width();

    double phi_factor = 1.0 + 0.1*sin(phi);

    return 2*rp/equil.bounds.get_r_width()*phi_factor/r;
}

double get_analytic_Bphi(const Equilibrium& equil, double r, double z, double phi){
    double phi_factor = 1.0 + 0.1*cos(phi);
    return phi_factor;
}


// Create a psi eq grid based off the analytic psi
Kokkos::View<double***, Kokkos::LayoutRight, HostType> analytic_3d_psi_grid(const Equilibrium& equil,
                                                                           const View<double*, HostType>& r,
                                                                           const View<double*, HostType>& z,
                                                                           const View<double*, HostType>& phi){

    Kokkos::View<double***, Kokkos::LayoutRight, HostType> psi_grid("psi_grid",phi.size(), z.size(), r.size());
    for (int k=0; k<phi.size(); k++){
        for (int j=0; j<z.size(); j++){
            for (int i=0; i<r.size(); i++){
                psi_grid(k,j,i) = get_circular_analytic_psi3d(equil, r(i), z(j), phi(k));
            }
        }
    }
    return psi_grid;
}

// Create a Br eq grid based off the analytic Br
Kokkos::View<double***, Kokkos::LayoutRight, HostType> analytic_3d_Br_grid(const Equilibrium& equil,
                                                                           const View<double*, HostType>& r,
                                                                           const View<double*, HostType>& z,
                                                                           const View<double*, HostType>& phi){

    Kokkos::View<double***, Kokkos::LayoutRight, HostType> Br_grid("Br_grid",phi.size(), z.size(), r.size());
    for (int k=0; k<phi.size(); k++){
        for (int j=0; j<z.size(); j++){
            for (int i=0; i<r.size(); i++){
                Br_grid(k,j,i) = get_analytic_Br(equil, r(i), z(j), phi(k));
            }
        }
    }
    return Br_grid;
}

// Create a Bz eq grid based off the analytic Bz
Kokkos::View<double***, Kokkos::LayoutRight, HostType> analytic_3d_Bz_grid(const Equilibrium& equil,
                                                                           const View<double*, HostType>& r,
                                                                           const View<double*, HostType>& z,
                                                                           const View<double*, HostType>& phi){

    Kokkos::View<double***, Kokkos::LayoutRight, HostType> Bz_grid("Bz_grid",phi.size(), z.size(), r.size());
    for (int k=0; k<phi.size(); k++){
        for (int j=0; j<z.size(); j++){
            for (int i=0; i<r.size(); i++){
                Bz_grid(k,j,i) = get_analytic_Bz(equil, r(i), z(j), phi(k));
            }
        }
    }
    return Bz_grid;
}

// Create a Bphi eq grid based off the analytic Bphi
Kokkos::View<double***, Kokkos::LayoutRight, HostType> analytic_3d_Bphi_grid(const Equilibrium& equil,
                                                                           const View<double*, HostType>& r,
                                                                           const View<double*, HostType>& z,
                                                                           const View<double*, HostType>& phi){

    Kokkos::View<double***, Kokkos::LayoutRight, HostType> Bphi_grid("Bphi_grid",phi.size(), z.size(), r.size());
    for (int k=0; k<phi.size(); k++){
        for (int j=0; j<z.size(); j++){
            for (int i=0; i<r.size(); i++){
                Bphi_grid(k,j,i) = get_analytic_Bphi(equil, r(i), z(j), phi(k));
            }
        }
    }
    return Bphi_grid;
}


}
#endif