magnetic_equil_files.hpp

Go to the documentation of this file.

#ifndef MAGNETIC_EQUIL_FILES_HPP
#define MAGNETIC_EQUIL_FILES_HPP

#include "space_settings.hpp"
#include "NamelistReader.hpp"
#include "globals.hpp"
#include "range_view.hpp"
#include "rz_bounds.hpp"
#include "file_reader.hpp"

extern "C" void xgc_read_eq_m3dc1_dims(int* mpsi, int* mr, int* mz, double* min_r, double* max_r, double* min_z, double* max_z);
extern "C" void xgc_read_eq_m3dc1(int mpsi, int mr, int mz, double min_r, double max_r, double min_z, double max_z,
                                  double* axis_r, double* axis_z, double* axis_b, double* x_psi, double* x_r, double* x_z,
                                  double* psi_grid, double* rgrid, double* zgrid, double* I, double* psi_rz);

struct MagneticEquilFiles{
    int mpsi;
    int mr;
    int mz;
    RZBounds bounds;
    double axis_r;
    double axis_z;
    double axis_b;
    double x_psi;
    double x_r;
    double x_z;
    View<double*, CLayout, HostType> psi_grid;
    View<double*, HostType> rgrid;
    View<double*, HostType> zgrid;
    View<double*, CLayout, HostType> I;
    View<double**, CLayout, HostType> psi_rz;

    // Reads and communicates eq dimensions
    void read_eq_dims(FileReader& file_reader){
        file_reader.skip_line(); // Skip first line (Header)
        file_reader.read(&mr, 1);
        file_reader.read(&mz, 1);
        file_reader.read(&mpsi, 1);
        file_reader.read(&bounds.min_r, 1);
        file_reader.read(&bounds.max_r, 1);
        file_reader.read(&bounds.min_z, 1);
        file_reader.read(&bounds.max_z, 1);
    }

    void read_eq(FileReader& file_reader){
        file_reader.read(&axis_r, 1);
        file_reader.read(&axis_z, 1);
        file_reader.read(&axis_b, 1);
        file_reader.read(&x_psi, 1);
        file_reader.read(&x_r, 1);
        file_reader.read(&x_z, 1);

        file_reader.read(psi_grid.data(), psi_grid.size());
        file_reader.read(I.data(), I.size());
        file_reader.read(psi_rz.data(), psi_rz.size());

        int end_flag;
        file_reader.read(&end_flag, 1);
        if(end_flag!=-1){
            exit_XGC("\nError reading eqd file: Did not encounter end flag (-1) where expected. The file format may be incorrect.\n");
        } 
    }

    // Read the files
    MagneticEquilFiles(bool read_m3dc1, const std::string& eq_filename){
        // file_reader if using simple eqd file
        FileReader file_reader;

        // First read in the dimensions
        if(is_rank_zero()){
            if(read_m3dc1){
                xgc_read_eq_m3dc1_dims(&mpsi, &mr, &mz, &bounds.min_r, &bounds.max_r, &bounds.min_z, &bounds.max_z);
            }else{
                file_reader.open(eq_filename);
                read_eq_dims(file_reader);
            }
        }

        // Broadcast dimensions to other processes
#ifdef USE_MPI
        int root_rank = 0;
        MPI_Bcast(&mr,   1, MPI_INT, root_rank, SML_COMM_WORLD);
        MPI_Bcast(&mz,   1, MPI_INT, root_rank, SML_COMM_WORLD);
        MPI_Bcast(&mpsi, 1, MPI_INT, root_rank, SML_COMM_WORLD);

        MPI_Bcast(&bounds, 4, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
#endif

        // Allocate views
        psi_grid = View<double*, CLayout, HostType>("psi_grid", mpsi);
        I = View<double*, CLayout, HostType>("I", mpsi);
        psi_rz = View<double**, CLayout, HostType>("psi_rz", mz, mr);

        // Set rgrid and zgrid to be linear
        rgrid = create_range_view("rgrid",bounds.min_r, bounds.max_r, mr);
        zgrid = create_range_view("zgrid",bounds.min_z, bounds.max_z, mz);

        // Read in data
        if(is_rank_zero()){
            if(read_m3dc1){
                xgc_read_eq_m3dc1(mpsi, mr, mz, bounds.min_r, bounds.max_r, bounds.min_z, bounds.max_z, &axis_r, &axis_z, &axis_b, &x_psi, &x_r, &x_z,
                            psi_grid.data(), rgrid.data(), zgrid.data(), I.data(), psi_rz.data());
            }else{
                read_eq(file_reader);
            }
        }

        // Broadcast values to other processes
#ifdef USE_MPI
        MPI_Bcast(&axis_r,1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast(&axis_z,1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast(&axis_b,1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);

        MPI_Bcast(&x_psi, 1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast(&x_r,   1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast(&x_z,   1, MPI_DOUBLE, root_rank, SML_COMM_WORLD);

        MPI_Bcast(psi_grid.data(), psi_grid.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast(I.data(),        I.size(),        MPI_DOUBLE, root_rank, SML_COMM_WORLD);
        MPI_Bcast(psi_rz.data(),   psi_rz.size(),   MPI_DOUBLE, root_rank, SML_COMM_WORLD);
#endif
    }

    // Analytic constructor
    MagneticEquilFiles(NLReader::NamelistReader& nlr){

        if(is_rank_zero()) printf("\nCreating a circular analytic magnetic equilibrium.\n");

        axis_r = 1.5;
        if(nlr.namelist_present("analytic_grid_param")){
            nlr.use_namelist("analytic_grid_param");
            axis_r = nlr.get<double>("axis_r", 1.5);
        }

        axis_z = 0.0;
        axis_b = 1.0;

        //nlr.use_namelist("analytic_eq_files_param");
        mpsi=100;
        mr=100;
        mz=100;
        bounds.min_r = axis_r - 1.0;
        bounds.max_r = axis_r + 1.0;
        bounds.min_z = axis_z - 1.0;
        bounds.max_z = axis_z + 1.0;
        x_psi = 2.0; // outside bounds
        x_r = axis_r;
        x_z = bounds.min_z - 1.0; // outside bounds

        // Set rgrid and zgrid to be linear
        rgrid = create_range_view("rgrid",bounds.min_r, bounds.max_r, mr);
        zgrid = create_range_view("zgrid",bounds.min_z, bounds.max_z, mz);
        psi_grid = create_range_view("psi_grid",0.0, 1.0, mpsi);

        // Allocate views
        I = View<double*, CLayout, HostType>("I", mpsi);
        Kokkos::deep_copy(I, 10000.0); // High safety factor for now - make input later

        // Set up psi_grid - concentric linear psi
        psi_rz = View<double**, CLayout, HostType>("psi_rz", mz, mr);
        for (int j=0; j<zgrid.size(); j++){
            for (int i=0; i<rgrid.size(); i++){
                double rp = (rgrid(i) - axis_r);
                double zp = (zgrid(j) - axis_z);
                psi_rz(j,i) = sqrt(rp*rp + zp*zp);
            }
        }
    }

    /* Returns MagneticEquilFiles object using the file-reader constructor. This may be updated so the input file
     * can specify an analytic grid.
     * */
    static MagneticEquilFiles eq_files_from_namelist(NLReader::NamelistReader& nlr){
        nlr.use_namelist("sml_param");

        if(nlr.get<bool>("sml_create_analytic_equil", false)==true) return MagneticEquilFiles(nlr);

        // Invert psi values if appropriate
        bool bt_sign = nlr.get<double>("sml_bt_sign", -1.0);
        bool bp_sign = nlr.get<double>("sml_bp_sign", 1.0);

        // Read file names and path from input file
        bool read_eqd = nlr.get<bool>("sml_read_eq", true);
        bool read_m3dc1 = nlr.get<bool>("sml_read_eq_m3dc1", false);
        // Ignore read_eqd if read_m3dc1 is true

        if(!(read_eqd || read_m3dc1)) exit_XGC("Error: must use either sml_read_eq or sml_read_eq_m3dc1");

        std::string input_file_dir = nlr.get<std::string>("sml_input_file_dir","./");

        nlr.use_namelist("eq_param");
        std::string eq_filename = nlr.get<std::string>("eq_filename","example_file.eqd"); 

        // Add path to file name
        eq_filename = input_file_dir + "/" + eq_filename;

        // Read file
        auto eq_file = MagneticEquilFiles(read_m3dc1, eq_filename);

        eq_file.x_psi *= bp_sign;
        eq_file.axis_b *= bt_sign;

        for(int i=0; i<eq_file.psi_grid.size(); i++) eq_file.psi_grid(i) *= bp_sign;
        for(int i=0; i<eq_file.I.size(); i++) eq_file.I(i) *= bt_sign;
        for(int i=0; i<eq_file.psi_rz.extent(0); i++)
            for(int j=0; j<eq_file.psi_rz.extent(1); j++)
                eq_file.psi_rz(i,j) *= bp_sign;

        // Return result
        return eq_file;
    }

};

#endif