Skip to content

Commit

Permalink
[TASK] non linear kicker integration
Browse files Browse the repository at this point in the history
work in progress
  • Loading branch information
PierreSchnizer committed Oct 10, 2023
1 parent 1d3a484 commit 44ed5d8
Show file tree
Hide file tree
Showing 11 changed files with 494 additions and 48 deletions.
4 changes: 2 additions & 2 deletions python/src/elements.cc
Original file line number Diff line number Diff line change
Expand Up @@ -410,8 +410,8 @@ void py_thor_scsi_init_elements(py::module &m)
TemplatedClasses<tsc::TpsaVariantType> templated_classes_tpsa
(m, "Tpsa");
auto classical_magnet_tpsa = templated_classes_tpsa.buildClasses(elem_type, field_kick_tpsa);
classical_magnet_add_set_methods<double, tse::ClassicalMagnetWithKnob<tsc::StandardDoubleType>> (classical_magnet_tpsa);
classical_magnet_add_set_methods<gtpsa::tpsa, tse::ClassicalMagnetWithKnob<tsc::TpsaVariantType>> (classical_magnet_tpsa);
classical_magnet_add_set_methods<double, tse::ClassicalMagnetWithKnob<tsc::StandardDoubleType>> (classical_magnet_tpsa);
classical_magnet_add_set_methods<gtpsa::tpsa, tse::ClassicalMagnetWithKnob<tsc::TpsaVariantType >> (classical_magnet_tpsa);



Expand Down
30 changes: 26 additions & 4 deletions python/src/interpolation.cc
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,13 @@
#include <gtpsa/python/utils.h>
#include <thor_scsi/core/field_interpolation.h>
#include <thor_scsi/core/multipoles.h>
#include <thor_scsi/custom/nonlinear_kicker_interpolation.h>
#include <complex>

namespace py = pybind11;
namespace gpy = gtpsa::python;
namespace tsc = thor_scsi::core;
namespace tsu = thor_scsi::custom;

typedef std::complex<double> cdbl;
typedef std::complex<double> cdbl_intern;
Expand Down Expand Up @@ -307,6 +309,21 @@ void py_thor_scsi_init_field_interpolation(py::module &m) {
/* air coil magnets helpers */
/* for communicating the field description */
PYBIND11_NUMPY_DTYPE(tsu::aircoil_filament_t, x, y, current);
py::class_<tsu::aircoil_filament_t> aircoil_filament(m, "aircoil_filament");
aircoil_filament
.def_readwrite("x", &tsu::aircoil_filament_t::x)
.def_readwrite("y", &tsu::aircoil_filament_t::y)
.def_readwrite("current", &tsu::aircoil_filament_t::y)
.def("__repr__", [](tsu::aircoil_filament_t& filament){
std::stringstream strm;
filament.show(strm, 10);
return strm.str();
})
.def(py::init([](const double x, const double y, const double current){
tsu::aircoil_filament_t f = {x, y, current};
return f;
}))
;

py::class_<
tsu:: AirCoilMagneticFieldKnobbed<tsc::StandardDoubleType>,
Expand All @@ -317,18 +334,23 @@ void py_thor_scsi_init_field_interpolation(py::module &m) {
aircoil_magnetic_field
.def("set_scale", &tsu::AirCoilMagneticFieldKnobbed<tsc::StandardDoubleType>::setScale)
.def("get_scale", &tsu::AirCoilMagneticFieldKnobbed<tsc::StandardDoubleType>::getScale)
.def("get_used_filaments", [](tsu::AirCoilMagneticFieldKnobbed<tsc::StandardDoubleType>& inst){
auto tmp = inst.getUsedFilaments();
return (tmp);
})
.def(py::init([](py::array_t<tsu::aircoil_filament_t, py::array::c_style|py::array::forcecast>& a) {
return tsu::AirCoilMagneticField(aircoil_filaments_from_array(a));
}) , "initalise with filaments (x, y, current)");

py::class_<
tsu:: NonlinearKickerKnobbed<tsc::StandardDoubleType>,
std::shared_ptr<tsu:: NonlinearKickerKnobbed<tsc::StandardDoubleType>>
tsu::NonLinearKickerInterpolationKnobbed<tsc::StandardDoubleType>,
std::shared_ptr<tsu::NonLinearKickerInterpolationKnobbed<tsc::StandardDoubleType>>
>
nonlinear_kicker(m, "NonlinearKicker", aircoil_magnetic_field);
nonlinear_kicker(m, "NonLinearKickerInterpolation", aircoil_magnetic_field);

nonlinear_kicker
.def(py::init([](py::array_t<tsu::aircoil_filament_t, py::array::c_style|py::array::forcecast>& a) {
return tsu::NonlinearKicker(aircoil_filaments_from_array(a));
return tsu::NonLinearKickerInterpolationKnobbed<tsc::StandardDoubleType>(aircoil_filaments_from_array(a));
}), "initalise with filaments (x, y, current) of one quater");

}
17 changes: 17 additions & 0 deletions src/thor_scsi/custom/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -31,3 +31,20 @@ target_link_libraries(test_nonlinear_kicker_interpolation
${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
)
add_test(nonlinear_kicker_interpolation test_nonlinear_kicker_interpolation)

add_executable(test_nonlinear_kicker test_nonlinear_kicker.cc)
target_include_directories(test_nonlinear_kicker
PUBLIC
"$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
)
target_link_libraries(test_nonlinear_kicker
thor_scsi
thor_scsi_core
gtpsa::gtpsa-c++
# PRIVATE Boost::boost
# quadmath
${Boost_PRG_EXEC_MONITOR_LIBRARY}
${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}
)
add_test(nonlinear_kicker test_nonlinear_kicker)
80 changes: 54 additions & 26 deletions src/thor_scsi/custom/aircoil_interpolation.h
Original file line number Diff line number Diff line change
@@ -1,9 +1,12 @@
#ifndef _THOR_SCSI_CUSTOM_AIRCOIL_INTERPOLATION_H_
#define _THOR_SCSI_CUSTOM_AIRCOIL_INTERPOLATION_H_ 1

#include <ostream>
#include <vector>
#include <gtpsa/utils.hpp>
#include <thor_scsi/core/field_interpolation.h>


namespace thor_scsi::custom {
struct aircoil_filament {
double x;
Expand All @@ -19,42 +22,67 @@ namespace thor_scsi::custom {
typedef struct aircoil_filament aircoil_filament_t;

template<class C>
class AirCoilMagneticFieldKnobbed : thor_scsi::core::Field2DInterpolationKnobbed<C> {
const std::vector<aircoil_filament_t> m_filaments;
class AirCoilMagneticFieldKnobbed : public thor_scsi::core::Field2DInterpolationKnobbed<C> {
const std::vector<aircoil_filament_t> m_filaments;
double m_scale;

protected:
// subclasses will "multiply the number of filaments during initalisation"
inline void setFilaments(const std::vector<aircoil_filament_t> filaments){ this->m_filaments = filaments;}

public:
AirCoilMagneticFieldKnobbed(const std::vector<aircoil_filament_t> filaments, const double scale=1e0)
AirCoilMagneticFieldKnobbed(const std::vector<aircoil_filament_t> filaments, const double scale=0e0)
: m_filaments(filaments)
, m_scale(scale)
{}

virtual ~AirCoilMagneticFieldKnobbed() {}

void setScale(const double scale) {this->m_scale = scale; }
const double getScale(void) const { return this->m_scale; }
void setScale(const double scale) { this->m_scale = scale; }
const double getScale(void) const { return this->m_scale; }

// required for inspection: e.g. if python transfers filaments correctly
inline std::vector<aircoil_filament_t> getUsedFilaments(void) const {
return this->m_filaments;
}


/*
* equivalent python code
r2 = np.sum(dpos ** 2, axis=0)
dx, dy = dpos
dpos_cross = np.array([dy, dx])
B = self.precomp * 1 / r2 * dpos_cross
B = np.sum(B, axis=1)
field[0] = B[0]
field[1] = B[1]
*/
* r2 = np.sum(dpos ** 2, axis=0)
* dx, dy = dpos
* dpos_cross = np.array([dy, dx])
* B = self.precomp * 1 / r2 * dpos_cross
* B = np.sum(B, axis=1)
* field[0] = B[0]
* field[1] = B[1]
*/
template<typename T>
inline void _field(const T& x, const T& y, T *Bx, T *By) const {
const double mu0 = 4 * M_PI * 1e-7;
const double precomp = mu0 / (2 * M_PI) * this->m_scale;
*Bx = 0e0;
*By = 0e0;
for(const auto& f: this->m_filaments){
const auto dx = x - f.x;
const auto dy = y - f.y;
const auto r2 = dx*dx + dy*dy;
*By += precomp * f.current / r2 * dx;
*Bx += precomp * f.current / r2 * dy;
}
inline void _field(const T& x, const T& y, T *pBx, T *pBy) const {
const double mu0 = 4 * M_PI * 1e-7;
const double precomp = mu0 / (2 * M_PI) * this->m_scale;

/*
* better allocate new ones and copy back
* otherwise: for tpsa objects don't forget to clear them
* including gradients!
*
* new allocation does that for us ...
*/
T Bx = gtpsa::same_as_instance(*pBx);
T By = gtpsa::same_as_instance(*pBx);

for(const auto& f: this->m_filaments){
const auto dx = x - f.x;
const auto dy = y - f.y;
const auto r2 = dx*dx + dy*dy;
By += precomp * f.current / r2 * dx;
Bx += precomp * f.current / r2 * dy;
}

*pBx = std::move(Bx);
*pBy = std::move(By);

}

template<typename T>
Expand All @@ -68,7 +96,7 @@ namespace thor_scsi::custom {
inline void field(const gtpsa::tpsa& x, const gtpsa::tpsa& y, gtpsa::tpsa *Bx, gtpsa::tpsa *By) const override final
{
this->_field(x, y, Bx, By);
}
}

inline void gradient(const double& x, const double& y, double *Gx, double *Gy) const override final
{ /* this->gradient(x, y, Gx, Gy); */ };
Expand Down
2 changes: 2 additions & 0 deletions src/thor_scsi/custom/nonlinear_kicker.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
// check that syntax is correct
#include <thor_scsi/custom/nonlinear_kicker.h>
83 changes: 83 additions & 0 deletions src/thor_scsi/custom/nonlinear_kicker.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
#ifndef _THOR_SCSI_ELEMENTS_NONLINEAR_KICKER_H_
#define _THOR_SCSI_ELEMENTS_NONLINEAR_KICKER_H_

#include <thor_scsi/elements/field_kick.h>
#include <thor_scsi/custom/nonlinear_kicker_interpolation.h>

#include <iomanip>

namespace thor_scsi::custom {
template<class C, typename = typename C::complex_type, typename = typename C::double_type> // CK for complex knob
class NonLinearKickerTypeWithKnob : public thor_scsi::elements::FieldKickKnobbed<C> {
protected:
using complex_type = typename C::complex_type;
using double_type = typename C::double_type;
using nonlinearkicker_knobbed = thor_scsi::custom::NonLinearKickerInterpolationKnobbed<C>;
public:
/**
*
* \verbatim embed:rst:leading-asterisk
*
* .. Warning::
*
* fix memory handling of interpolation object
* add move constructor
*
* \endverbatim
*/
inline NonLinearKickerTypeWithKnob(const Config &config) : thor_scsi::elements::FieldKickKnobbed<C>(config){
std::vector<thor_scsi::custom::aircoil_filament_t> off {{200e0, 200e0, 0e0}};
auto tmp = std::make_shared<nonlinearkicker_knobbed>(off);
this->intp = std::move(std::dynamic_pointer_cast<thor_scsi::core::Field2DInterpolationKnobbed<C>>(tmp));
/*
std::cerr << "NonLinearKicker Type: " << std::setw(10) << name << " interpolation object " << this->intp << std::endl;
auto parent = dynamic_cast<FieldKick*>(this);
std::cerr << "FieldKick: interpolation object " << parent->getFieldInterpolator() << std::endl;
*/
}

const char* type_name(void) const override { return "NonLinearKicker"; };

inline auto getFieldInterpolator(void) const {
std::stringstream strm;
strm << __FILE__ << ":mpole.getFieldInterpolator: this->intp " << static_cast<void *>(this->intp.get());
auto p = std::dynamic_pointer_cast<nonlinearkicker_knobbed>(this->intp);
strm << " p " << static_cast<void *>(p.get());
if(!p){
strm<<" cast failed!";
throw std::runtime_error(strm.str());
}
auto cp = std::const_pointer_cast<nonlinearkicker_knobbed>(p);
strm << " cp " << static_cast<void *>(cp.get());
if(!cp){
strm<<" cast failed!";
throw std::runtime_error(strm.str());
}
// std::cerr << strm.str() << std::endl;
return cp;
}

inline void setFieldInterpolator(std::shared_ptr<thor_scsi::custom::NonLinearKickerInterpolationKnobbed<C>> intp) {
auto p = std::dynamic_pointer_cast<thor_scsi::core::Field2DInterpolationKnobbed<C>>(intp);
if(!p){
throw std::runtime_error("Could not cast multipole to Field2DInterpolation");
}
this->intp = p;
}

//thor_scsi::core::TwoDimensionalMultipoles* intp;
};


typedef class NonLinearKickerTypeWithKnob<core::StandardDoubleType> NonLinearKickerType;
typedef class NonLinearKickerTypeWithKnob<core::TpsaVariantType> NonLinearKickerTypeTpsa;

} // Name space

#endif // _THOR_SCSI_ELEMENTS_NONLINEAR_KICKER_H_
/*
* Local Variables:
* mode: c++
* c-file-style: "python"
* End:
*/
8 changes: 4 additions & 4 deletions src/thor_scsi/custom/nonlinear_kicker_interpolation.cc
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,12 @@ const std::vector<tsu::aircoil_filament_t> tsu::construct_aircoil_filaments(


template<class C>
void tsu::NonlinearKickerKnobbed<C>::show(std::ostream& strm, int level) const
void tsu::NonLinearKickerInterpolationKnobbed<C>::show(std::ostream& strm, int level) const
{
strm << "NonlinearKickerKnobbed({ ";
strm << "NonLinearKickerInterpolationKnobbed({ ";
tsu::AirCoilMagneticFieldKnobbed<C>::show(strm, level);
strm << "})";
}

template void tsu::NonlinearKickerKnobbed<tsc::StandardDoubleType>::show(std::ostream& strm, int level) const;
template void tsu::NonlinearKickerKnobbed<tsc::TpsaVariantType>::show(std::ostream& strm, int level) const;
template void tsu::NonLinearKickerInterpolationKnobbed<tsc::StandardDoubleType>::show(std::ostream& strm, int level) const;
template void tsu::NonLinearKickerInterpolationKnobbed<tsc::TpsaVariantType>::show(std::ostream& strm, int level) const;
9 changes: 6 additions & 3 deletions src/thor_scsi/custom/nonlinear_kicker_interpolation.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,17 +15,20 @@ namespace thor_scsi::custom {
const std::vector<aircoil_filament_t> construct_aircoil_filaments(const std::vector<aircoil_filament_t>& filaments_one_quater);

template<class C>
class NonlinearKickerKnobbed: public AirCoilMagneticFieldKnobbed<C> {
class NonLinearKickerInterpolationKnobbed: public AirCoilMagneticFieldKnobbed<C> {

public:
NonlinearKickerKnobbed(const std::vector<aircoil_filament_t>& filaments_one_quater)
NonLinearKickerInterpolationKnobbed(const std::vector<aircoil_filament_t> filaments_one_quater)
: AirCoilMagneticFieldKnobbed<C>(construct_aircoil_filaments(filaments_one_quater))
{}

inline void setFilaments(const std::vector<aircoil_filament_t> filaments_one_quater) {
this->setFilaments(construct_aircoil_filaments(filaments_one_quater));
}
void show(std::ostream&, int level) const override;
};

typedef NonlinearKickerKnobbed<thor_scsi::core::StandardDoubleType> NonlinearKicker;
typedef NonLinearKickerInterpolationKnobbed<thor_scsi::core::StandardDoubleType> NonLinearKickerInterpolation;
} //namespace thor_scsi::custom

#endif /* _THOR_SCSI_CUSTOM_NONLINEAR_KICKER_INTERPOLATION_H_ */
Loading

0 comments on commit 44ed5d8

Please sign in to comment.