Commit 0f48558c authored by cianciosa's avatar cianciosa Committed by Cianciosa, Mark
Browse files

Add implementation for mpex fields.

parent 513bea6a

graph_driver/mpex_example.sh

0 → 100755
+1 −0
Original line number Diff line number Diff line 
./graph_driver/xrays --absorption_model=weak_damping --dispersion=ordinary_wave --endtime=0.15 --equilibrium=mpex --equilibrium_file=../graph_tests/mpex.nc --mpex_ps1_current=4500 --mpex_ps2_current=4500 --mpex_tr1_current=600 --mpex_tr2_current=160 --mpex_ne_scale=1.0E19 --mpex_te_scale=1000 --init_kx --init_kx_mean=-700.0 --init_ky_dist=normal --init_ky_mean=-100.0 --init_ky_sigma=10.0 --init_kz_dist=normal --init_kz_mean=0.0 --init_kz_sigma=10.0 --init_w_dist=normal --init_w_mean=700 --init_w_sigma=10.0 --init_x_mean=0.2 --init_y_dist=normal --init_y_mean=0.0 --init_y_sigma=0.05 --init_z_dist=normal --init_z_mean=2.0 --init_z_sigma=0.0 --num_rays=1 --num_times=100000 --solver=rk4 --sub_steps=100 --use_cyl_xy

graph_driver/xrays.cpp

+22 −1
Original line number Diff line number Diff line
@@ -378,6 +378,14 @@ equilibrium::shared<T, SAFE_MATH> make_equilibrium(const commandline::parser &cl 


    if (eq == "efit") {

        return equilibrium::make_efit<T, SAFE_MATH> (file_name);

    } else if (eq == "mpex") {

        return equilibrium::make_mpex<T, SAFE_MATH> (file_name,

                                                     cl.get_option_value<T> ("mpex_ne_scale"),

                                                     cl.get_option_value<T> ("mpex_te_scale"),

                                                     cl.get_option_value<T> ("mpex_ps1_current"),

                                                     cl.get_option_value<T> ("mpex_ps2_current"),

                                                     cl.get_option_value<T> ("mpex_tr1_current"),

                                                     cl.get_option_value<T> ("mpex_tr2_current"));

    } else {

        return equilibrium::make_vmec<T, SAFE_MATH> (file_name);

    }
@@ -976,13 +984,21 @@ commandline::parser parse_commandline(int argc, const char * argv[]) { 
    });

    cl.add_option("equilibrium",       true,  "Equilibrium to use.", {

        "efit",

        "vmec"

        "vmec",

        "mpex"

    });

    cl.add_option("equilibrium_file",  true,  "File to read the equilibrium from.");

    cl.add_option("init_w_dist",       true,  "Distribution function for wave frequency.", {

        "uniform",

        "normal"

    });

    cl.add_option("mpex_ne_scale",     true,  "Scale factor for electron density profiles.");

    cl.add_option("mpex_te_scale",     true,  "Scale factor for electron temperature profiles.");

    cl.add_option("mpex_ps2_current",  true,  "Current for the ps2 coil set.");

    cl.add_option("mpex_ps1_current",  true,  "Current for the ps1 coil set.");

    cl.add_option("mpex_ps2_current",  true,  "Current for the ps2 coil set.");

    cl.add_option("mpex_tr1_current",  true,  "Current for the tr1 coil set.");

    cl.add_option("mpex_tr2_current",  true,  "Current for the tr2 coil set.");

    cl.add_option("init_w_mean",       true,  "Mean value for the wave frequency distribution function.");

    cl.add_option("init_w_sigma",      true,  "Standard deviation of for the wave frequency distribution function.");

    cl.add_option("init_kx_dist",      true,  "Distribution function for wave number in the x direction.", {
@@ -1097,6 +1113,11 @@ int main(int argc, const char * argv[]) { 


    typedef double base;



    std::cout << "Using " << cl.get_option_value<std::string> ("equilibrium") << " equilibrium from " << cl.get_option_value<std::string> ("equilibrium_file") << std::endl;

    std::cout << "Using " << cl.get_option_value<std::string> ("dispersion") << " dispersion relation" << std::endl;

    std::cout << "Using " << cl.get_option_value<std::string> ("solver") <<  " solver methd" << std::endl;

    std::cout << "Using " << cl.get_option_value<std::string> ("absorption_model") << " absorption model" << std::endl << std::endl;



    trace_ray<base> (cl, num_times, sub_steps, num_rays);

    calculate_power<std::complex<base>, use_safe_math> (cl,

                                                        num_times,

graph_framework/equilibrium.hpp

+747 −157

File changed.

Preview size limit exceeded, changes collapsed.

graph_framework/output.hpp

+9 −0
Original line number Diff line number Diff line
@@ -23,7 +23,14 @@ namespace output { 
///  @param[in] status Error status code.

//------------------------------------------------------------------------------

    static void check_error(const int status) {

#ifdef NDEBUG

        if (status) {

            std::cerr << nc_strerror(status) << std::endl;

            exit(status);

        }

#else

        assert(status == NC_NOERR && nc_strerror(status));

#endif

    }



//------------------------------------------------------------------------------
@@ -85,7 +92,9 @@ namespace output { 
///  @brief Destructor.

//------------------------------------------------------------------------------

        ~result_file() {

            sync.lock();

            check_error(nc_close(ncid));

            sync.unlock();

        }



//------------------------------------------------------------------------------

graph_korc/xkorc.cpp

+29 −14
Original line number Diff line number Diff line 
#include <random>



#include "../graph_framework/equilibrium.hpp"

#include "../graph_framework/timing.hpp"

#include "../graph_framework/output.hpp"
@@ -11,7 +13,7 @@ template<jit::float_scalar T> 
void run_korc() {

    const timing::measure_diagnostic t_total("Total Time");

    

    const size_t num_particles = 1000000;

    const size_t num_particles = 100000;

    std::cout << "Num particles " << num_particles << std::endl;

    std::vector<std::thread> threads(std::max(std::min(static_cast<unsigned int> (jit::context<T>::max_concurrency()),

                                                       static_cast<unsigned int> (num_particles)),
@@ -40,22 +42,34 @@ void run_korc() { 
            

            std::cout << "Local num particles " << local_num_particles << std::endl;



            std::mt19937_64 engine((thread_number + 1)*static_cast<uint64_t> (std::chrono::system_clock::to_time_t(std::chrono::system_clock::now())));

            std::normal_distribution<T> normal(static_cast<T> (0.0), static_cast<T> (0.1));



            auto ux = graph::variable<T> (local_num_particles, "u_{x}");

            auto uy = graph::variable<T> (local_num_particles, "u_{y}");

            auto uz = graph::variable<T> (local_num_particles, "u_{z}");



            ux->set(0.0);

            uy->set(0.99);

            uz->set(0.1);

            uy->set(0.9);

            for (size_t j = 0; j < local_num_particles; j++) {

                T temp = normal(engine);

                uz->set(j, temp);

                while (0.9*0.9 + temp*temp > 1.0) {

                    temp = normal(engine);

                    uz->set(j, temp);

                }

            }

            

            auto x = graph::variable<T> (local_num_particles, "x");

            auto y = graph::variable<T> (local_num_particles, "y");

            auto z = graph::variable<T> (local_num_particles, "z");

            auto pos = graph::vector(x, y, z);

            

            x->set(1.7);

            y->set(0.0);

            z->set(0.0);

            for (size_t j = 0; j < local_num_particles; j++) {

                x->set(j, 1.7 + normal(engine));

                y->set(j, 0.0);

                z->set(j, normal(engine));

            }

            

            auto u_vec = graph::vector(ux, uy, uz);
@@ -143,19 +157,20 @@ void run_korc() { 
            

            const timing::measure_diagnostic t_run("Run Time");

            work.pre_run();

            for (size_t i = 0; i < 1000000; i++) {

/*                sync.join();

            for (size_t j = 0; j < 10000; j++) {

                sync.join();

                work.wait();

                sync = std::thread([&file, &dataset] () -> void {

                    dataset.write(file);

                });*/

                

                });

                for (size_t k = 0; k < 10; k++) {

                    work.run();

                }

            }

            work.wait();



            sync.join();

            dataset.write(file);

            //dataset.write(file);

            work.wait();

            

            t_run.print();