diff --git a/graph_framework/cuda_context.hpp b/graph_framework/cuda_context.hpp index 2f371fb04334454d667cc55e4662d1c985cd5621..044b4cb954eb6f17337094215de6f96e0d08544b 100644 --- a/graph_framework/cuda_context.hpp +++ b/graph_framework/cuda_context.hpp @@ -746,6 +746,7 @@ namespace gpu { graph::output_nodes &outputs, graph::map_nodes &setters, jit::register_map ®isters, + jit::register_map &indices, const jit::register_usage &usage) { for (auto &[out, in] : setters) { graph::shared_leaf a = out->compile(source_buffer, diff --git a/graph_korc/xkorc.cpp b/graph_korc/xkorc.cpp index 95ddfdbd825b6a53cb4945efa9fc09b53dbf87b5..30a9e24a61bb41a7530085f7d4dd99db738ad2b0 100644 --- a/graph_korc/xkorc.cpp +++ b/graph_korc/xkorc.cpp @@ -2,79 +2,75 @@ #include "../graph_framework/timing.hpp" //------------------------------------------------------------------------------ -/// @brief Main program of the driver. +/// @brief Run Korc /// -/// @param[in] argc Number of commandline arguments. -/// @param[in] argv Array of commandline arguments. +/// @tparam T Base type. //------------------------------------------------------------------------------ -int main(int argc, const char * argv[]) { - START_GPU - (void)argc; - (void)argv; - +template +void run_korc() { const timeing::measure_diagnostic t_total("Total Time"); - - const size_t num_particles = 10000000; + + const size_t num_particles = 1; std::cout << "Num particles " << num_particles << std::endl; - std::vector threads(std::max(std::min(static_cast (jit::context::max_concurrency()), + std::vector threads(std::max(std::min(static_cast (jit::context::max_concurrency()), static_cast (num_particles)), static_cast (1))); - + const size_t batch = num_particles/threads.size(); const size_t extra = num_particles%threads.size(); - + for (size_t i = 0, ie = threads.size(); i < ie; i++) { threads[i] = std::thread([num_particles, batch, extra] (const size_t thread_number) -> void { const size_t local_num_particles = batch + (extra > thread_number ? 1 : 0); - + const timeing::measure_diagnostic t_setup("Setup Time"); - - auto eq = equilibrium::make_efit (EFIT_FILE); - //auto eq = equilibrium::make_slab_density (); + + auto eq = equilibrium::make_efit (EFIT_FILE); + //auto eq = equilibrium::make_slab_density (); auto b0 = eq->get_characteristic_field(thread_number); - const double q = 1.602176634E-19; - const double me = 9.1093837139E-31; - const double c = 299792458.0; - + const T q = 1.602176634E-19; + const T me = 9.1093837139E-31; + const T c = 299792458.0; + auto gryo_period = me/(q*b0); std::cout << "gryo_period " << gryo_period->evaluate().at(0) << std::endl; auto larmor_radius = c*gryo_period; std::cout << "larmor_radius " << larmor_radius->evaluate().at(0) << std::endl; - + std::cout << "Local num particles " << local_num_particles << std::endl; - - auto ux = graph::variable (local_num_particles, "u_{x}"); - auto uy = graph::variable (local_num_particles, "u_{y}"); - auto uz = graph::variable (local_num_particles, "u_{z}"); - - ux->set(0.99); - uy->set(0.0); - uz->set(0.0); - - auto x = graph::variable (local_num_particles, "x"); - auto y = graph::variable (local_num_particles, "y"); - auto z = graph::variable (local_num_particles, "z"); + + auto ux = graph::variable (local_num_particles, "u_{x}"); + auto uy = graph::variable (local_num_particles, "u_{y}"); + auto uz = graph::variable (local_num_particles, "u_{z}"); + + ux->set(0.0); + uy->set(0.99); + uz->set(0.1); + + auto x = graph::variable (local_num_particles, "x"); + auto y = graph::variable (local_num_particles, "y"); + auto z = graph::variable (local_num_particles, "z"); auto pos = graph::vector(x, y, z); - + x->set(1.7); y->set(0.0); z->set(0.0); - + auto u_vec = graph::vector(ux, uy, uz); - - auto gamma = graph::variable (local_num_particles, "\\gamma"); - - auto dt = graph::constant (0.25); - + + auto gamma = graph::variable (local_num_particles, "\\gamma"); + + auto dt = graph::constant (0.25); + auto gamma_init = 1.0/graph::sqrt(1.0 - u_vec->dot(u_vec)); - + auto u_init = gamma_init*u_vec; - + auto b_vec = eq->get_magnetic_field(pos->get_x(), pos->get_y(), pos->get_z())/b0; - - workflow::manager work(thread_number); + + workflow::manager work(thread_number); work.add_preitem({ graph::variable_cast(ux), graph::variable_cast(uy), @@ -86,25 +82,25 @@ int main(int argc, const char * argv[]) { {u_init->get_z(), graph::variable_cast(uz)}, {gamma_init, graph::variable_cast(gamma)} }, "initalize_gamma"); - + auto u_prime = u_vec - dt*u_vec->cross(b_vec)/(2.0*gamma); - + auto tau = -0.5*dt*b_vec; auto tau_sq = tau->dot(tau); auto speed_sq = u_prime->dot(u_prime); auto sigma = 1.0 + speed_sq - tau_sq; auto ustar = u_prime->dot(tau); - + auto gamma_next = graph::sqrt(0.5*(sigma + graph::sqrt(sigma*sigma + 4.0*(tau_sq + ustar*ustar)))); auto t = tau/gamma_next; - + auto s = 1.0 + t->dot(t); auto u_prime_dot_t = u_prime->dot(t); - + auto u_next = (u_prime + u_prime_dot_t*t + u_prime->cross(t))/s; - - auto pos_next = pos + larmor_radius*dt*u_next/gamma; - + + auto pos_next = pos + larmor_radius*dt*u_next/gamma_next; + work.add_item({ graph::variable_cast(x), graph::variable_cast(y), @@ -122,26 +118,69 @@ int main(int argc, const char * argv[]) { {u_next->get_z(), graph::variable_cast(uz)}, {gamma_next, graph::variable_cast(gamma)} }, "step"); - + work.compile(); + + std::ostringstream stream; + stream << "korc_" << thread_number << ".nc"; + + output::result_file file(stream.str(), local_num_particles); + output::data_set dataset(file); + + dataset.create_variable(file, "x", x, work.get_context()); + dataset.create_variable(file, "y", y, work.get_context()); + dataset.create_variable(file, "z", z, work.get_context()); + dataset.create_variable(file, "ux", ux, work.get_context()); + dataset.create_variable(file, "uy", uy, work.get_context()); + dataset.create_variable(file, "uz", uz, work.get_context()); + dataset.create_variable(file, "gamma", gamma, work.get_context()); + + file.end_define_mode(); + std::thread sync([]{}); + t_setup.print(); - + const timeing::measure_diagnostic t_run("Run Time"); work.pre_run(); for (size_t i = 0; i < 1000000; i++) { + /*sync.join(); + work.wait(); + sync = std::thread([&file, &dataset] () -> void { + dataset.write(file); + });*/ + work.run(); } work.wait(); + + sync.join(); + dataset.write(file); + work.wait(); + t_run.print(); }, i); } - + for (std::thread &t : threads) { t.join(); } - + std::cout << std::endl << "Timing:" << std::endl; t_total.print(); +} + +//------------------------------------------------------------------------------ +/// @brief Main program of the driver. +/// +/// @param[in] argc Number of commandline arguments. +/// @param[in] argv Array of commandline arguments. +//------------------------------------------------------------------------------ +int main(int argc, const char * argv[]) { + START_GPU + (void)argc; + (void)argv; + + run_korc (); END_GPU } diff --git a/graph_playground/xplayground.cpp b/graph_playground/xplayground.cpp index 2d5702f51e79bb573ca1411efd7919e3b35710fd..6f3855dc88e6ca10f2e102cea7daaec56dc73e9f 100644 --- a/graph_playground/xplayground.cpp +++ b/graph_playground/xplayground.cpp @@ -4,6 +4,10 @@ //------------------------------------------------------------------------------ #include "../graph_framework/jit.hpp" +#include + +#include "../graph_framework/equilibrium.hpp" + //------------------------------------------------------------------------------ /// @brief Main program of the playground. /// @@ -17,7 +21,47 @@ int main(int argc, const char * argv[]) { // Insert code here. No code should be commited to this file beyond this // template. + auto r = graph::variable (1000, "r"); + auto phi = graph::variable (1000, "\\phi"); + auto z = graph::variable (1000, "z"); + + auto eq = equilibrium::make_efit (EFIT_FILE); + + auto bvec = eq->get_magnetic_field(r*graph::cos(phi), + r*graph::sin(phi), + z); + + bvec->get_x()->to_latex(); + std::cout << "\\\\" << std::endl; + bvec->get_y()->to_latex(); + std::cout << "\\\\" << std::endl; + bvec->get_z()->to_latex(); + std::cout << "\\\\" << std::endl; + std::cout << std::endl << std::endl << std::endl; + + r->set(static_cast (1.7)); + z->set(static_cast (0.0)); + std::vector temp(1000); + for (size_t i = 0; i < 1000; i++) { + temp[i] = 2.0*std::numbers::pi_v*i/999.0; + } + phi->set(temp); + workflow::manager work(0); + work.add_item({ + graph::variable_cast(r), + graph::variable_cast(phi), + graph::variable_cast(z) + }, { + bvec->get_x(), + bvec->get_y(), + bvec->get_z() + }, {}, "bvec_kernel"); + work.compile(); + work.run(); + for (size_t i = 0; i < 1000; i++) { + work.print(i, {r, phi, z, bvec->get_x(), bvec->get_y(), bvec->get_z()}); + } END_GPU }