Commit e467dab8 authored by Cianciosa, Mark's avatar Cianciosa, Mark
Browse files

Allow the arguments of piecewise constants to be deferred. This fixed a bug... 

Allow the arguments of piecewise constants to be deferred. This fixed a bug where the EFIT test case was not reflecting off the O-Mode cutoff since the spline arguments were not correct.
parent 7f9ea865

graph_framework/CMakeLists.txt

+0 −2
Original line number Diff line number Diff line
@@ -41,8 +41,6 @@ target_precompile_headers (rays 
                           $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/cpu_context.hpp>

                           $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/piecewise.hpp>

                           $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/newton.hpp>

                           $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/coordinates.hpp>

                           $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/cublic_splines.hpp>

                           $<$<BOOL:${USE_METAL}>:$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/metal_context.hpp>>

                           $<$<BOOL:${USE_CUDA}>:$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/cuda_context.hpp>>

)

graph_framework/backend.hpp

+6 −0
Original line number Diff line number Diff line
@@ -28,6 +28,12 @@ namespace backend { 


    public:

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

///  @brief Construct an empty buffer backend.

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

        buffer() :

        memory() {}



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

///  @brief Construct a buffer backend with a size.

///

///  @params[in] s Size of he data buffer.

graph_framework/coordinates.hpp

deleted100644 → 0
+0 −79
Original line number Diff line number Diff line 
//------------------------------------------------------------------------------

///  @file coordinates.hpp

///  @brief Define graph generation routines for different coodinates.

///

///  Defines graphs for different coordinate systems.

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



#ifndef coordinates_h

#define coordinates_h



#include <array>



#include "vector.hpp"



namespace coordinates {

//******************************************************************************

//  Cartesian

//******************************************************************************

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

///  @brief Class representing cartesian coodinates.

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

    template<typename T>

    class cartesian {

    protected:

///  Base coordinates.

        std::array<graph::shared_leaf<T>, 3> base;

        

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

///  @brief Construct cartesian coordinates.

///

///  @params[in] size    Number of elements in the variable.

///  @params[in] symbols Symbols for the coordinates.

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

        cartesian(const size_t size,

                  const std::array<std::string, 3> symbols) :

        base({graph::variable<T> (size, symbols.at(0)),

              graph::variable<T> (size, symbols.at(1)),

              graph::variable<T> (size, symbols.at(2))}) {}

        

    public:

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

///  @brief Construct cartesian coordinates.

///

///  @params[in] size Number of elements in the variable.

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

        cartesian(const size_t size) :

        base({graph::variable<T> (size, "x"),

              graph::variable<T> (size, "y"),

              graph::variable<T> (size, "z")}) {}



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

///  @brief Convert to cartesian.

///

///  @returns The base in cartesian coordinates.

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

        virtual std::array<graph::shared_leaf<T>, 3> to_cartesian() const {

            return base;

        }



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

///  @brief Get basis vectors.

///

///  @returns The basis vectors in cartesian coordinates.

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

        virtual std::array<graph::shared_vector<T>, 3> get_basis() {

            return {graph::vector(graph::one<T> (),

                                  graph::one<T> (),

                                  graph::one<T> ()),

                    graph::vector(graph::one<T> (),

                                  graph::one<T> (),

                                  graph::one<T> ()),

                    graph::vector(graph::one<T> (),

                                  graph::one<T> (),

                                  graph::one<T> ())};

        }

    };

}



#endif /* coordinates_h */

graph_framework/cpu_context.hpp

+19 −19
Original line number Diff line number Diff line
@@ -142,7 +142,7 @@ namespace gpu { 
                exit(1);

            }



            std::vector<T *> buffers;

            std::map<std::string, T *> buffers;



            for (auto &input : inputs) {

                if (!kernel_arguments.contains(input.get())) {
@@ -151,20 +151,20 @@ namespace gpu { 
                    memcpy(arg.data(), buffer.data(), buffer.size()*sizeof(T));

                    kernel_arguments[input.get()] = arg;

                }

                buffers.push_back(kernel_arguments[input.get()].data());

                buffers[jit::to_string('v', input.get())] = kernel_arguments[input.get()].data();

            }

            for (auto &output : outputs) {

                if (!kernel_arguments.contains(output.get())) {

                    std::vector<T> arg(num_rays);

                    kernel_arguments[output.get()] = arg;

                }

                buffers.push_back(kernel_arguments[output.get()].data());

                buffers[jit::to_string('o', output.get())] = kernel_arguments[output.get()].data();

            }



            std::cout << "  Function pointer: " << reinterpret_cast<size_t> (kernel) << std::endl;



            return [kernel, buffers] () mutable {

                ((void (*)(std::vector<T *> &))kernel)(buffers);

                ((void (*)(std::map<std::string, T *> &))kernel)(buffers);

            };

        }
@@ -241,7 +241,8 @@ namespace gpu { 
///  @params[in,out] source_buffer Source buffer stream.

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

        void create_header(std::stringstream &source_buffer) {

            source_buffer << "#include <vector>" << std::endl;

            source_buffer << "#include <map>" << std::endl;

            source_buffer << "#include <string>" << std::endl;

            if (jit::is_complex<T> ()) {

                source_buffer << "#include <complex>" << std::endl;

            } else {
@@ -269,18 +270,19 @@ namespace gpu { 
            source_buffer << std::endl;

            source_buffer << "extern \"C\" void " << name << "(" << std::endl;

            

            source_buffer << "    vector<";

            source_buffer << "    std::map<std::string, ";

            jit::add_type<T> (source_buffer);

            source_buffer << " *> &args) {" << std::endl;

            

            source_buffer << "    for (size_t i = 0; i < " << size << "; i++) {" << std::endl;

            for (size_t i = 0, ie = inputs.size(); i < ie; i++) {

                registers[inputs[i].get()] = jit::to_string('r', inputs[i].get());



            for (auto &input : inputs) {

                registers[input.get()] = jit::to_string('r', input.get());

                source_buffer << "        const ";

                jit::add_type<T> (source_buffer);

                source_buffer << " " << registers[inputs[i].get()]

                              << " = args[" << i << "][i]; //" << inputs[i]->get_symbol() << std::endl;

                arg_index[inputs[i].get()] = i;

                source_buffer << " " << registers[input.get()]

                              << " = args[std::string(\"" << jit::to_string('v', input.get())

                              << "\")][i]; //" << input->get_symbol() << std::endl;

            }

        }
@@ -298,15 +300,13 @@ namespace gpu { 
                                   jit::register_map &registers) {

            for (auto &[out, in] : setters) {

                graph::shared_leaf<T> a = out->compile(source_buffer, registers);

                source_buffer << "        args[" << arg_index[in.get()];

                source_buffer << "][i] = " << registers[a.get()] << ";" << std::endl;

                source_buffer << "        args[std::string(\"" << jit::to_string('v', in.get());

                source_buffer << "\")][i] = " << registers[a.get()] << ";" << std::endl;

            }



            for (size_t i = 0, ie = outputs.size(); i < ie; i++) {

                graph::shared_leaf<T> a = outputs[i]->compile(source_buffer, registers);

                source_buffer << "        args[" << arg_index.size() + i

                              << "][i] = " << registers[a.get()] << ";"

                              << std::endl;

            for (auto &out : outputs) {

                graph::shared_leaf<T> a = out->compile(source_buffer, registers);

                source_buffer << "        args[std::string(\"" << jit::to_string('o', out.get());

                source_buffer << "\")][i] = " << registers[a.get()] << ";" << std::endl;

            }

                    

            source_buffer << "    }" << std::endl;

graph_framework/cublic_splines.hpp

deleted100644 → 0
+0 −136
Original line number Diff line number Diff line 
//------------------------------------------------------------------------------

///  @file cublic\_splines.hpp

///  @brief Creates function for cubic spline interpolated quantities.

///

///  Defines a tree of operations that allows automatic differentiation.

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



#ifndef cublic_splines_h

#define cublic_splines_h



#include <vector>



#include "arithmetic.hpp"



namespace spline {

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

///  @brief Build a 1D spline.

///

///  @params[in] min_x     Minimum argument.

///  @params[in] dx        Grid spacing.

///  @params[in] c0_buffer Zeroth coefficient buffer.

///  @params[in] c1_buffer First coefficient buffer.

///  @params[in] c2_buffer Second coefficient buffer.

///  @params[in] c3_buffer Third coefficient buffer.

///  @params[in] x         Non-normalized function argument.

///  @returns The graph of a 1D spline function.

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

    template<typename T>

    graph::shared_leaf<T> make_1D(graph::shared_leaf<T> min_x,

                                  graph::shared_leaf<T> dx,

                                  const std::vector<T> c0_buffer,

                                  const std::vector<T> c1_buffer,

                                  const std::vector<T> c2_buffer,

                                  const std::vector<T> c3_buffer,

                                  graph::shared_leaf<T> x) {

        auto arg_norm = (x - min_x)/dx;



        auto c0 = graph::piecewise_1D(c0_buffer, arg_norm);

        auto c1 = graph::piecewise_1D(c1_buffer, arg_norm);

        auto c2 = graph::piecewise_1D(c2_buffer, arg_norm);

        auto c3 = graph::piecewise_1D(c3_buffer, arg_norm);

        

        return c0 + c1*arg_norm + c2*arg_norm*arg_norm + c3*arg_norm*arg_norm*arg_norm;

    }



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

///  @brief Build a 2D spline.

///

///  @params[in] min_x      Minimum argument.

///  @params[in] dx         Grid spacing.

///  @params[in] min_y      Minimum argument.

///  @params[in] dy         Grid spacing.

///  @params[in] c00_buffer Zeroth coefficient buffer.

///  @params[in] c01_buffer First coefficient buffer.

///  @params[in] c02_buffer Second coefficient buffer.

///  @params[in] c03_buffer Third coefficient buffer.

///  @params[in] c10_buffer Zeroth coefficient buffer.

///  @params[in] c11_buffer First coefficient buffer.

///  @params[in] c12_buffer Second coefficient buffer.

///  @params[in] c13_buffer Third coefficient buffer.

///  @params[in] c20_buffer Zeroth coefficient buffer.

///  @params[in] c21_buffer First coefficient buffer.

///  @params[in] c22_buffer Second coefficient buffer.

///  @params[in] c23_buffer Third coefficient buffer.

///  @params[in] c30_buffer Zeroth coefficient buffer.

///  @params[in] c31_buffer First coefficient buffer.

///  @params[in] c32_buffer Second coefficient buffer.

///  @params[in] c33_buffer Third coefficient buffer.

///  @params[in] x          Non-normalized function argument.

///  @params[in] y          Non-normalized function argument.

///  @returns The graph of a 1D spline function.

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

    template<typename T>

    graph::shared_leaf<T> make_2D(graph::shared_leaf<T> min_x,

                                  graph::shared_leaf<T> dx,

                                  graph::shared_leaf<T> min_y,

                                  graph::shared_leaf<T> dy,

                                  const std::vector<T> c00_buffer,

                                  const std::vector<T> c01_buffer,

                                  const std::vector<T> c02_buffer,

                                  const std::vector<T> c03_buffer,

                                  const std::vector<T> c10_buffer,

                                  const std::vector<T> c11_buffer,

                                  const std::vector<T> c12_buffer,

                                  const std::vector<T> c13_buffer,

                                  const std::vector<T> c20_buffer,

                                  const std::vector<T> c21_buffer,

                                  const std::vector<T> c22_buffer,

                                  const std::vector<T> c23_buffer,

                                  const std::vector<T> c30_buffer,

                                  const std::vector<T> c31_buffer,

                                  const std::vector<T> c32_buffer,

                                  const std::vector<T> c33_buffer,

                                  graph::shared_leaf<T> x,

                                  graph::shared_leaf<T> y,

                                  const size_t num_columns) {

        auto x_norm = (x - min_x)/dx;

        auto y_norm = (y - min_y)/dy;



        auto c00 = graph::piecewise_2D(c00_buffer, num_columns, x_norm, y_norm);

        auto c01 = graph::piecewise_2D(c01_buffer, num_columns, x_norm, y_norm);

        auto c02 = graph::piecewise_2D(c02_buffer, num_columns, x_norm, y_norm);

        auto c03 = graph::piecewise_2D(c03_buffer, num_columns, x_norm, y_norm);

        auto c10 = graph::piecewise_2D(c10_buffer, num_columns, x_norm, y_norm);

        auto c11 = graph::piecewise_2D(c11_buffer, num_columns, x_norm, y_norm);

        auto c12 = graph::piecewise_2D(c12_buffer, num_columns, x_norm, y_norm);

        auto c13 = graph::piecewise_2D(c13_buffer, num_columns, x_norm, y_norm);

        auto c20 = graph::piecewise_2D(c20_buffer, num_columns, x_norm, y_norm);

        auto c21 = graph::piecewise_2D(c21_buffer, num_columns, x_norm, y_norm);

        auto c22 = graph::piecewise_2D(c22_buffer, num_columns, x_norm, y_norm);

        auto c23 = graph::piecewise_2D(c23_buffer, num_columns, x_norm, y_norm);

        auto c30 = graph::piecewise_2D(c30_buffer, num_columns, x_norm, y_norm);

        auto c31 = graph::piecewise_2D(c31_buffer, num_columns, x_norm, y_norm);

        auto c32 = graph::piecewise_2D(c32_buffer, num_columns, x_norm, y_norm);

        auto c33 = graph::piecewise_2D(c33_buffer, num_columns, x_norm, y_norm);

        

        return c00 +

               c01*y_norm +

               c02*y_norm*y_norm +

               c03*y_norm*y_norm*y_norm +

               c10*x_norm +

               c11*x_norm*y_norm +

               c12*x_norm*y_norm*y_norm +

               c13*x_norm*y_norm*y_norm*y_norm +

               c20*x_norm*x_norm +

               c21*x_norm*x_norm*y_norm +

               c22*x_norm*x_norm*y_norm*y_norm +

               c23*x_norm*x_norm*y_norm*y_norm*y_norm +

               c30*x_norm*x_norm*x_norm +

               c31*x_norm*x_norm*x_norm*y_norm +

               c32*x_norm*x_norm*x_norm*y_norm*y_norm +

               c33*x_norm*x_norm*x_norm*y_norm*y_norm*y_norm;

    }

}



#endif /* cublic_splines_h */