Commit 63a58c38 authored by cianciosa's avatar cianciosa
Browse files

Commit changes before merging in bug fix.

parent 36bd4ae9

graph_framework/vector.hpp

+18 −0
Original line number Diff line number Diff line
@@ -294,6 +294,24 @@ namespace graph { 
                      l->get_z() + r->get_z());

    }



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

///  @brief Subtraction operator.

///

///  @tparam T         Base type of the calculation.

///  @tparam SAFE_MATH Use safe math operations.

///

///  @param[in] l Left vector.

///  @param[in] r Right vector.

///  @returns The vector vector addition.

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

    template<jit::float_scalar T, bool SAFE_MATH=false>

    shared_vector<T, SAFE_MATH> operator-(shared_vector<T, SAFE_MATH> l,

                                          shared_vector<T, SAFE_MATH> r) {

        return vector(l->get_x() - r->get_x(),

                      l->get_y() - r->get_y(),

                      l->get_z() - r->get_z());

    }



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

///  @brief Multiplication operator.

///

graph_korc/CMakeLists.txt

+1 −1
Original line number Diff line number Diff line 
add_tool_target (xkorc)



if (${USE_PCH})

    target_precompile_headers (xrays_bench REUSE_FROM xrays)

    target_precompile_headers (xkorc REUSE_FROM xrays)

endif ()

graph_korc/xkorc.cpp

+50 −36
Original line number Diff line number Diff line
@@ -14,7 +14,7 @@ int main(int argc, const char * argv[]) { 


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



    const size_t num_particles = 1000000;

    const size_t num_particles = 1;

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

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

                                                       static_cast<unsigned int> (num_particles)),
@@ -64,13 +64,17 @@ int main(int argc, const char * argv[]) { 


            auto gamma = graph::variable<double> (local_num_particles, "\\gamma");



            auto dt = graph::constant<double> (0.1);

            auto dt = graph::constant<double> (0.01);



            auto gamma_init = graph::sqrt(1.0 - ux*ux - uy*uy - uz*uz);



            auto u_init = gamma_init*u_vec;



            workflow::manager<double> work(0);

            auto b_vec = eq->get_magnetic_field(pos->get_x(),

                                                pos->get_y(),

                                                pos->get_z())/b0;



            workflow::manager<double> work(thread_number);

            work.add_preitem({

                graph::variable_cast(x),

                graph::variable_cast(y),
@@ -86,27 +90,23 @@ int main(int argc, const char * argv[]) { 
                {gamma_init, graph::variable_cast(gamma)}

            }, "initalize_gamma");



            auto pos_next = pos + larmor_radius*dt*u_vec/gamma;

            

            auto b_vec = eq->get_magnetic_field(pos_next->get_x(),

                                                pos_next->get_y(),

                                                pos_next->get_z())/b0;

            

            auto u_prime = u_vec + dt*u_vec->cross(b_vec)/(2.0*gamma);

            auto u_prime = u_vec - dt*u_vec->cross(b_vec)/(2.0*gamma);



            auto tau = dt*0.5*b_vec;

            auto tau = -0.5*dt*b_vec;

            auto tau_sq = tau->dot(tau);

            auto speed_sq = u_vec->dot(u_vec);

            auto speed_sq = u_prime->dot(u_prime);

            auto sigma = 1.0 + speed_sq - tau_sq;

            auto ustar = u_vec->dot(tau);

            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/(1.0 + t->dot(t));

            auto s = 1.0 + t->dot(t);

            auto u_prime_dot_t = u_prime->dot(t);



            auto u_next = s*(u_prime + u_prime_dot_t*t + u_prime->cross(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;



            work.add_item({

                graph::variable_cast(x),
@@ -116,7 +116,11 @@ int main(int argc, const char * argv[]) { 
                graph::variable_cast(uy),

                graph::variable_cast(uz),

                graph::variable_cast(gamma)

            }, {}, {

            }, {

                tau_sq,

                tau->get_x()*tau->get_x(),

                tau->get_y()*tau->get_y(),

            }, {

                {pos_next->get_x(), graph::variable_cast(x)},

                {pos_next->get_y(), graph::variable_cast(y)},

                {pos_next->get_z(), graph::variable_cast(z)},
@@ -125,14 +129,24 @@ int main(int argc, const char * argv[]) { 
                {u_next->get_z(), graph::variable_cast(uz)},

                {gamma_next, graph::variable_cast(gamma)}

            }, "step");

            tau->get_x()->to_latex();

            std::cout << "\\\\" << std::endl;

            tau->get_y()->to_latex();

            std::cout << "\\\\" << std::endl;

            (tau->get_x()*tau->get_x())->to_latex();

            std::cout << "\\\\" << std::endl;

            (tau->get_y()*tau->get_y())->to_latex();

            std::cout << "\\\\" << std::endl;

    

            work.compile();

            t_setup.print();



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

            work.pre_run();

            work.print(0, {x, y, z, ux, uy, uz, gamma, tau_sq, tau->get_x()*tau->get_x(), tau->get_y()*tau->get_y()});

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

                work.run();

                work.print(0, {x, y, z, ux, uy, uz, gamma, tau_sq, tau->get_x()*tau->get_x(), tau->get_y()*tau->get_y()});

            }

            work.wait();

            t_run.print();