Merge pull request #13 from martinschwinzerl/kparasch_master

Kparasch master

Settings.cmake.default

@@ -39,8 +39,7 @@ option( SIXTRACKL_ENABLE_BENCHMARK_TESTS "Enable benchmark tests"         OFF )
 
 set( SIXTRACKL_TRACK_BEAMBEAM4D             "enabled" CACHE STRING "" )
 set( SIXTRACKL_TRACK_BEAMBEAM6D             "enabled" CACHE STRING "" )
-set( SIXTRACKL_TRACK_SPACECHARGE_COASTING   "enabled" CACHE STRING "" )
-set( SIXTRACKL_TRACK_SPACECHARGE_BUNCHED    "enabled" CACHE STRING "" )
+set( SIXTRACKL_TRACK_SPACECHARGE            "enabled" CACHE STRING "" )
 set( SIXTRACKL_TRACK_TRICUB                 "skip"    CACHE STRING "" )
 
 # Aperture check features: allowed values are always, never, conditional

cmake/SetupFeatures.cmake

@@ -6,8 +6,7 @@ if(  NOT SETUP_SIXTRL_FEATURES_FINISHED )
 
     set_property( CACHE SIXTRACKL_TRACK_BEAMBEAM4D
                         SIXTRACKL_TRACK_BEAMBEAM6D
-                        SIXTRACKL_TRACK_SPACECHARGE_COASTING
-                        SIXTRACKL_TRACK_SPACECHARGE_BUNCHED
+                        SIXTRACKL_TRACK_SPACECHARGE
                         SIXTRACKL_TRACK_TRICUB
                   PROPERTY STRINGS enabled disabled skip )
 
@@ -47,37 +46,19 @@ if(  NOT SETUP_SIXTRL_FEATURES_FINISHED )
         endif()
     endif()
 
-    set_property( CACHE SIXTRACKL_TRACK_SPACECHARGE_COASTING
-                  PROPERTY HELPSTRING
-                  "Track over coasting space-charge beam element" )
-
-    if( SIXTRACKL_TRACK_SPACECHARGE_COASTING )
-        if( "${SIXTRACKL_TRACK_SPACECHARGE_COASTING}" STREQUAL "skip" )
-            set( SIXTRL_TRACK_SC_COASTING_FLAG "${SIXTRL_TRACK_MAP_SKIP_VALUE}" )
-            set( SIXTRL_TRACK_SC_COASTING_FLAG_STR "skip" )
-        elseif( "${SIXTRACKL_TRACK_SPACECHARGE_COASTING}" STREQUAL "disabled" )
-            set( SIXTRL_TRACK_SC_COASTING_FLAG "${SIXTRL_TRACK_MAP_DISABLED_VALUE}" )
-            set( SIXTRL_TRACK_SC_COASTING_FLAG_STR "disabled" )
-        else()
-            set( SIXTRL_TRACK_SC_COASTING_FLAG "${SIXTRL_TRACK_MAP_ENABLED_VALUE}" )
-            set( SIXTRL_TRACK_SC_COASTING_FLAG_STR "enabled" )
-        endif()
-    endif()
+    set_property( CACHE SIXTRACKL_TRACK_SPACECHARGE PROPERTY HELPSTRING
+                  "Track over frozen space-charge beam element" )
 
-    set_property( CACHE SIXTRACKL_TRACK_SPACECHARGE_BUNCHED
-                  PROPERTY HELPSTRING
-                  "Track over bunched space-charge beam element" )
-
-    if( SIXTRACKL_TRACK_SPACECHARGE_BUNCHED )
-        if( "${SIXTRACKL_TRACK_SPACECHARGE_BUNCHED}" STREQUAL "skip" )
-            set( SIXTRL_TRACK_SC_BUNCHED_FLAG "${SIXTRL_TRACK_MAP_SKIP_VALUE}" )
-            set( SIXTRL_TRACK_SC_BUNCHED_FLAG_STR "skip" )
-        elseif( "${SIXTRACKL_TRACK_SPACECHARGE_BUNCHED}" STREQUAL "disabled" )
-            set( SIXTRL_TRACK_SC_BUNCHED_FLAG "${SIXTRL_TRACK_MAP_DISABLED_VALUE}" )
-            set( SIXTRL_TRACK_SC_BUNCHED_FLAG_STR "disabled" )
+    if( SIXTRACKL_TRACK_SPACECHARGE )
+        if( "${SIXTRACKL_TRACK_SPACECHARGE}" STREQUAL "skip" )
+            set( SIXTRL_TRACK_SC_FLAG "${SIXTRL_TRACK_MAP_SKIP_VALUE}" )
+            set( SIXTRL_TRACK_SC_FLAG_STR "skip" )
+        elseif( "${SIXTRACKL_TRACK_SPACECHARGE}" STREQUAL "disabled" )
+            set( SIXTRL_TRACK_SC_FLAG "${SIXTRL_TRACK_MAP_DISABLED_VALUE}" )
+            set( SIXTRL_TRACK_SC_FLAG_STR "disabled" )
         else()
-            set( SIXTRL_TRACK_SC_BUNCHED_FLAG "${SIXTRL_TRACK_MAP_ENABLED_VALUE}" )
-            set( SIXTRL_TRACK_SC_BUNCHED_FLAG_STR "enabled" )
+            set( SIXTRL_TRACK_SC_FLAG "${SIXTRL_TRACK_MAP_ENABLED_VALUE}" )
+            set( SIXTRL_TRACK_SC_FLAG_STR "enabled" )
         endif()
     endif()
 
@@ -105,10 +86,8 @@ if(  NOT SETUP_SIXTRL_FEATURES_FINISHED )
           set( SIXTRL_TRACK_BEAMBEAM4D_FLAG_STR  \"${SIXTRL_TRACK_BEAMBEAM4D_FLAG_STR}\" )
           set( SIXTRL_TRACK_BEAMBEAM6D_FLAG      \"${SIXTRL_TRACK_BEAMBEAM4D_FLAG}\" )
           set( SIXTRL_TRACK_BEAMBEAM6D_FLAG_STR  \"${SIXTRL_TRACK_BEAMBEAM4D_FLAG_STR}\" )
-          set( SIXTRL_TRACK_SC_COASTING_FLAG     \"${SIXTRL_TRACK_SC_COASTING_FLAG}\" )
-          set( SIXTRL_TRACK_SC_COASTING_FLAG_STR \"${SIXTRL_TRACK_SC_COASTING_FLAG_STR}\" )
-          set( SIXTRL_TRACK_SC_BUNCHED_FLAG      \"${SIXTRL_TRACK_SC_BUNCHED_FLAG}\" )
-          set( SIXTRL_TRACK_SC_BUNCHED_FLAG_STR  \"${SIXTRL_TRACK_SC_BUNCHED_FLAG_STR}\" )
+          set( SIXTRL_TRACK_SPACECHARGE_FLAG     \"${SIXTRL_TRACK_SPACECHARGE_FLAG}\" )
+          set( SIXTRL_TRACK_SPACECHARGE_FLAG_STR \"${SIXTRL_TRACK_SPACECHARGE_FLAG_STR}\" )
           set( SIXTRL_TRACK_TRICUB_FLAG          \"${SIXTRL_TRACK_TRICUB_FLAG}\"
           set( SIXTRL_TRACK_TRICUB_FLAG_STR      \"${SIXTRL_TRACK_TRICUB_FLAG_STR}\" )" )
 

python/sixtracklib/beam_elements.py

@@ -95,15 +95,130 @@ def set_ksl(self, value, order):
 
 
 class RFMultipole(CObject):
-    _typeid = 128
+    _typeid = 256  # This is subject to change
     order = CField(0, 'int64', default=0, const=True, alignment=8)
     voltage = CField(1, 'real', default=0.0, alignment=8)
     frequency = CField(2, 'real', default=0.0, alignment=8)
     lag = CField(3, 'real', default=0.0, alignment=8)
-    knl = CField(4, 'real', pointer=True, length='order+1', alignment=8)
-    ksl = CField(5, 'real', pointer=True, length='order+1', alignment=8)
-    pn = CField(6, 'real', pointer=True, length='order+1', alignment=8),
-    ps = CField(7, 'real', pointer=True, length='order+1', alignment=8),
+    bal = CField(4, 'real', pointer=True, length='2*order+2',
+                    default=0.0, alignment=8)
+    phase = CField(5, 'real', pointer=True, length='2*order+2',
+                    default=0.0, alignment=8)
+
+    def __init__(
+            self,
+            order=None,
+            knl=None,
+            ksl=None,
+            pn=None,
+            ps=None,
+            bal=None,
+            p=None,
+            **kwargs):
+        if bal is None and \
+                (knl is not None or ksl is not None or
+                 pn is not None or ps is not None or order is not None):
+            if knl is None:
+                knl = []
+            if ksl is None:
+                ksl = []
+            if pn is None:
+                pn = []
+            if ps is None:
+                ps = []
+            if order is None:
+                order = 0
+
+            n = max((order + 1), max(len(knl), len(ksl), len(pn), len(ps)))
+            assert(n > 0)
+
+            _knl = np.array(knl)
+            nknl = np.zeros(n, dtype=_knl.dtype)
+            nknl[:len(knl)] = knl
+            knl = nknl
+            del(_knl)
+            assert(len(knl) == n)
+
+            _ksl = np.array(ksl)
+            nksl = np.zeros(n, dtype=_ksl.dtype)
+            nksl[:len(ksl)] = ksl
+            ksl = nksl
+            del(_ksl)
+            assert(len(ksl) == n)
+
+            _pn = np.array(pn)
+            npn = np.zeros(n, dtype=_pn.dtype)
+            npn[:len(pn)] = pn
+            pn = npn
+            del(_pn)
+            assert(len(pn) == n)
+
+            _ps = np.array(ps)
+            nps = np.zeros(n, dtype=_ps.dtype)
+            nps[:len(ps)] = ps
+            ps = nps
+            del(_ps)
+            assert(len(ps) == n)
+
+            order = n - 1
+            bal = np.zeros(2 * order + 2)
+            p = np.zeros(2 * order + 2)
+
+            idx = np.array([ii for ii in range(0, len(knl))])
+            inv_factorial = 1.0 / factorial(idx, exact=True)
+            bal[0::2] = knl * inv_factorial
+            bal[1::2] = ksl * inv_factorial
+
+            p[0::2] = pn
+            p[1::2] = ps
+
+            kwargs["bal"] = bal
+            kwargs["phase"] = p
+            kwargs["order"] = order
+
+        elif bal is not None and bal and len(bal) > 2 and ((len(bal) % 2) == 0)\
+                and p is not None and p and len(p) > 2 and ((len(p) % 2) == 0):
+            kwargs["bal"] = bal
+            kwargs["phase"] = p
+            kwargs["order"] = (len(bal) - 2) / 2
+
+        super().__init__(**kwargs)
+
+    @property
+    def knl(self):
+        idx = np.array([ii for ii in range(0, len(self.bal), 2)])
+        return self.bal[idx] * factorial(idx // 2, exact=True)
+
+    @property
+    def ksl(self):
+        idx = np.array([ii for ii in range(0, len(self.bal), 2)])
+        return self.bal[idx + 1] * factorial(idx // 2, exact=True)
+
+    def set_knl(self, value, order):
+        assert order <= self.order
+        self.bal[order * 2] = value / factorial(order, exact=True)
+
+    def set_ksl(self, value, order):
+        assert order <= self.order
+        self.bal[order * 2 + 1] = value / factorial(order, exact=True)
+
+    @property
+    def pn(self):
+        idx = np.array([ii for ii in range(0, len(self.p), 2)])
+        return self.phase[idx]
+
+    @property
+    def ps(self):
+        idx = np.array([ii for ii in range(0, len(self.p), 2)])
+        return self.phase[idx + 1]
+
+    def set_pn(self, value, order):
+        assert order <= self.order
+        self.phase[order * 2] = value
+
+    def set_ps(self, value, order):
+        assert order <= self.order
+        self.phase[order * 2 + 1] = value
 
 
 class Cavity(CObject):
@@ -291,7 +406,6 @@ def set_half_axes_squ(self, a_squ, b_squ):
         self.a_b_squ = a_squ * b_squ
         return self
 
-
 class LimitZeta(CObject):
     _typeid = 13
     min_zeta = CField(0, 'float64', default=-1e18, alignment=8)
@@ -443,6 +557,59 @@ def __init__(self, **kwargs):
         else:
             super().__init__(**kwargs)
 
+class LimitRectEllipse(CObject):
+    _typeid = 16
+    max_x = CField(0, 'float64', default=+1.0, alignment=8)
+    max_y = CField(1, 'float64', default=+1.0, alignment=8)
+    a_squ = CField(2, 'float64', default=+1.0, alignment=8)
+    b_squ = CField(3, 'float64', default=+1.0, alignment=8)
+    a_b_squ = CField(4, 'float64', alignment=8)
+
+    def __init__(
+            self,
+            max_x=None,
+            max_y=None,
+            a_squ=None,
+            b_squ=None,
+            **kwargs):
+        if max_x is None:
+            max_x = 1.0
+        if max_y is None:
+            max_y = 1.0
+        if a_squ is None and 'a' in kwargs:
+            a = kwargs.get('a')
+            if a is not None and a > 0.0:
+                a_squ = a * a
+        if a_squ is None:
+            a_squ = 1.0
+
+        if b_squ is None and 'b' in kwargs:
+            b = kwargs.get('b')
+            if b is not None and b > 0.0:
+                b_squ = b * b
+        if b_squ is None:
+            b_squ = 1.0
+
+        if max_x < 0.0:
+            raise ValueError("max_x has to be positive definite")
+
+        if max_y < 0.0:
+            raise ValueError("max_y has to be_positive definite")
+
+        if a_squ < 0.0 or b_squ < 0.0:
+            raise ValueError("a_squ and b_squ have to be positive definite")
+
+        super().__init__(max_x=max_x, max_y=max_y, a_squ=a_squ, b_squ=b_squ,
+                         a_b_squ=a_squ * b_squ, **kwargs)
+
+    def set_half_axes(self, a, b):
+        return self.set_half_axes_squ(a * a, b * b)
+
+    def set_half_axes_squ(self, a_squ, b_squ):
+        self.a_squ = a_squ
+        self.b_squ = b_squ
+        self.a_b_squ = a_squ * b_squ
+        return self
 
 class DipoleEdge(CObject):
     _typeid = 64
@@ -500,13 +667,19 @@ class Elements(object):
                      'BeamBeam6D': BeamBeam6D,
                      'SpaceChargeCoasting': SpaceChargeCoasting,
                      'SpaceChargeBunched': SpaceChargeBunched,
+                     'LimitRect': LimitRect,
+                     'LimitEllipse': LimitEllipse,
+                     'LimitRectEllipse': LimitRectEllipse,
                      'DipoleEdge': DipoleEdge,
                      #                     'Line': Line,
-                     # 'TriCub': TriCub,
                      # 'Solenoid': Solenoid,
                      # 'Wire': Wire,
                      }
 
+    @staticmethod
+    def add_element_type(cls, cls_name):
+        Elements.element_types.update({cls_name: cls})
+
     def _mk_fun(self, buff, cls):
         def fun(*args, **nargs):
             # print(cls.__name__,nargs)

python/sixtracklib/config.py.template

@@ -23,8 +23,7 @@ SHARED_LIB = os.path.join( os.path.dirname( os.path.realpath( __file__ ) ),
 
 TRACK_BEAMBEAM4D  = "@SIXTRL_TRACK_BEAMBEAM4D_FLAG_STR@"
 TRACK_BEAMBEAM6D  = "@SIXTRL_TRACK_BEAMBEAM6D_FLAG_STR@"
-TRACK_SC_COASTING = "@SIXTRL_TRACK_SC_COASTING_FLAG_STR@"
-TRACK_SC_BUNCHED  = "@SIXTRL_TRACK_SC_BUNCHED_FLAG_STR@"
+TRACK_SPACECHARGE = "@SIXTRL_TRACK_SC_FLAG_STR@"
 TRACK_TRICUB      = "@SIXTRL_TRACK_TRICUB_FLAG_STR@"
 
 # -----------------------------------------------------------------------------

python/sixtracklib/control.py

@@ -25,38 +25,40 @@
     st_Controller_get_remap_cobject_buffer_debug_kernel_id, \
     st_Controller_set_remap_cobject_buffer_debug_kernel_id, \
     st_Controller_get_num_available_nodes, \
-    st_Controller_has_selected_node, st_Controller_get_default_node_info_base, \
-    st_Controller_get_default_node_index, st_Controller_get_default_node_id, \
-    st_Controller_has_default_node, \
-    st_Controller_get_max_available_node_index, \
-    st_Controller_get_min_available_node_index, \
-    st_Controller_get_ptr_node_info_base_by_index, \
+    st_Controller_get_available_node_indices, \
+    st_Controller_is_node_available_by_index, \
     st_Controller_get_ptr_node_id_by_index, \
+    st_Controller_get_ptr_node_info_base_by_index, \
+    st_Controller_get_min_available_node_index, \
+    st_Controller_get_max_available_node_index, \
+    st_Controller_has_default_node, st_Controller_get_default_node_id, \
+    st_Controller_get_default_node_index, \
+    st_Controller_get_default_node_info_base, \
+    st_Controller_has_selected_node, \
     st_Controller_get_ptr_selected_node_id, \
-    st_Controller_get_selected_node_index, st_Controller_is_default_node, \
+    st_Controller_get_selected_node_index, \
     st_Controller_get_ptr_selected_node_info_base, \
-    st_Controller_can_change_selected_node, st_Controller_is_default_node_id, \
+    st_Controller_can_change_selected_node, \
     st_Controller_can_directly_change_selected_node, \
-    st_Controller_can_unselect_node, st_Controller_is_node_available_by_index, \
-    st_Controller_get_available_node_indices, \
-    st_Controller_is_default_node_index, \
+    st_Controller_can_unselect_node, st_Controller_is_default_node_id, \
+    st_Controller_is_default_node, st_Controller_is_default_node_index, \
+    st_Controller_is_default_platform_id_and_device_id, \
     st_Controller_get_ptr_node_info_base, \
-    st_Controller_get_ptr_node_info_base_by_index, \
     st_Controller_get_ptr_node_info_base_by_node_id, \
     st_Controller_get_ptr_node_info_base_by_platform_id_and_device_id, \
-    st_Controller_get_node_index, st_Controller_get_node_index_by_node_id, \
+    st_Controller_get_node_index_by_node_id, st_Controller_get_node_index, \
     st_Controller_get_node_index_by_platform_id_and_device_id, \
-    st_Controller_is_default_platform_id_and_device_id, \
-    st_Controller_is_node_available, st_Controller_is_node_available_by_index, \
     st_Controller_is_node_available_by_node_id, \
+    st_Controller_is_node_available, \
     st_Controller_is_node_available_by_platform_id_and_device_id, \
-    st_Controller_select_node, \
     st_Controller_select_node_by_index, st_Controller_select_node_by_node_id, \
+    st_Controller_select_node, \
     st_Controller_select_node_by_plaform_id_and_device_id, \
-    st_Controller_unselect_node, st_Controller_unselect_node_by_index, \
+    st_Controller_change_selected_node, st_Controller_unselect_node, \
+    st_Controller_unselect_node_by_index, \
     st_Controller_unselect_node_by_node_id, \
-    st_Controller_unselect_node_by_platform_id_and_device_id, \
-    st_Controller_change_selected_node
+    st_Controller_unselect_node_by_platform_id_and_device_id
+
 from .stcommon import st_Argument_delete, \
     st_Argument_get_arch_id, st_Argument_has_arch_string, \
     st_Argument_get_arch_string, st_Argument_send_again, \