Skip to content

Commit

Permalink
Merge pull request #129 from martinschwinzerl/hotfix/issue128_fix_tra…
Browse files Browse the repository at this point in the history 
…ck_line

Hotfix/issue128 fix track line
  • Loading branch information
@martinschwinzerl
martinschwinzerl authored Apr 12, 2020
2 parents b1e3fab + 6b78c92 commit 947643d
Show file tree
Hide file tree
Showing 6 changed files with 418 additions and 8 deletions.
10 changes: 5 additions & 5 deletions sixtracklib/common/be_limit/track.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,18 +72,18 @@ SIXTRL_INLINE NS(track_status_t) NS(Track_particle_limit_global)(

#endif /* SIXTRL_APERTURE_X_LIMIT && SIXTRL_APERTURE_Y_LIMIT */

index_t state = NS(Particles_get_state_value)( p, idx );

real_t const x = NS(Particles_get_x_value)( p, idx );
real_t const y = NS(Particles_get_y_value)( p, idx );

real_t const sign_x = ( real_t )( ( ZERO < x ) - ( real_t )( x < ZERO ) );
real_t const sign_y = ( real_t )( ( ZERO < y ) - ( real_t )( y < ZERO ) );

index_t const new_state = ( index_t )(
( ( sign_x * x ) < X_LIMIT ) & ( ( sign_y * y ) < Y_LIMIT ) );

SIXTRL_ASSERT( NS(Particles_is_not_lost_value)( p, idx ) );
NS(Particles_set_state_value)( p, idx, new_state );
state &= ( index_t )( ( ( sign_x * x ) < X_LIMIT ) &
( ( sign_y * y ) < Y_LIMIT ) );

NS(Particles_set_state_value)( p, idx, state );
return SIXTRL_TRACK_SUCCESS;
}

Expand Down
94 changes: 94 additions & 0 deletions sixtracklib/common/be_multipole/be_multipole.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,8 +146,35 @@ SIXTRL_FN SIXTRL_STATIC int NS(MultiPole_compare_values_with_treshold)(
SIXTRL_BE_ARGPTR_DEC const NS(MultiPole) *const SIXTRL_RESTRICT rhs,
NS(multipole_real_t) const treshold );

SIXTRL_FN SIXTRL_STATIC SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole) const*
NS(BufferIndex_get_const_multipole)(
SIXTRL_BUFFER_OBJ_ARGPTR_DEC const NS(Object) *const index_obj );

SIXTRL_FN SIXTRL_STATIC SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
NS(BufferIndex_get_multipole)( SIXTRL_BUFFER_OBJ_ARGPTR_DEC NS(Object)* index_obj );

SIXTRL_FN SIXTRL_STATIC SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole) const*
NS(BeamElements_managed_buffer_get_const_multipole)(
SIXTRL_BUFFER_DATAPTR_DEC unsigned char const* SIXTRL_RESTRICT pbuffer,
NS(buffer_size_t) const be_index, NS(buffer_size_t) const slot_size );

SIXTRL_FN SIXTRL_STATIC SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
NS(BeamElements_managed_buffer_get_multipole)(
SIXTRL_BUFFER_DATAPTR_DEC unsigned char* SIXTRL_RESTRICT pbuffer,
NS(buffer_size_t) const be_index, NS(buffer_size_t) const slot_size );

#if !defined( _GPUCODE )

SIXTRL_STATIC SIXTRL_HOST_FN SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole) const*
NS(BeamElements_buffer_get_const_multipole)(
SIXTRL_BUFFER_ARGPTR_DEC const NS(Buffer) *const SIXTRL_RESTRICT buffer,
NS(buffer_size_t) const be_index );

SIXTRL_STATIC SIXTRL_HOST_FN SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
NS(BeamElements_buffer_get_multipole)(
SIXTRL_BUFFER_ARGPTR_DEC NS(Buffer)* SIXTRL_RESTRICT buffer,
NS(buffer_size_t) const be_index );

SIXTRL_FN SIXTRL_STATIC bool NS(MultiPole_can_be_added)(
SIXTRL_BUFFER_ARGPTR_DEC const NS(Buffer) *const SIXTRL_RESTRICT buffer,
NS(multipole_order_t) const order,
Expand Down Expand Up @@ -841,8 +868,75 @@ SIXTRL_INLINE int NS(MultiPole_compare_values_with_treshold)(
return compare_value;
}

SIXTRL_INLINE SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole) const*
NS(BufferIndex_get_const_multipole)(
SIXTRL_BUFFER_OBJ_ARGPTR_DEC const NS(Object) *const index_obj )
{
typedef NS(MultiPole) beam_element_t;
typedef SIXTRL_BUFFER_OBJ_DATAPTR_DEC beam_element_t const* ptr_to_be_t;
ptr_to_be_t ptr_to_be = SIXTRL_NULLPTR;

if( ( index_obj != SIXTRL_NULLPTR ) &&
( NS(Object_get_type_id)( index_obj ) == NS(OBJECT_TYPE_MULTIPOLE) ) &&
( NS(Object_get_size)( index_obj ) >= sizeof( beam_element_t ) ) )
{
ptr_to_be = ( ptr_to_be_t )( uintptr_t
)NS(Object_get_begin_addr)( index_obj );
}

return ptr_to_be;
}

SIXTRL_INLINE SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
NS(BufferIndex_get_multipole)(
SIXTRL_BUFFER_OBJ_ARGPTR_DEC NS(Object)* index_obj )
{
return ( SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
)NS(BufferIndex_get_const_multipole)( index_obj );
}

SIXTRL_INLINE SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole) const*
NS(BeamElements_managed_buffer_get_const_multipole)(
SIXTRL_BUFFER_DATAPTR_DEC unsigned char const* SIXTRL_RESTRICT pbuffer,
NS(buffer_size_t) const be_index, NS(buffer_size_t) const slot_size )
{
return NS(BufferIndex_get_const_multipole)(
NS(ManagedBuffer_get_const_object)( pbuffer, be_index, slot_size ) );
}

SIXTRL_INLINE SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
NS(BeamElements_managed_buffer_get_multipole)(
SIXTRL_BUFFER_DATAPTR_DEC unsigned char* SIXTRL_RESTRICT pbuffer,
NS(buffer_size_t) const be_index, NS(buffer_size_t) const slot_size )
{
return NS(BufferIndex_get_multipole)(
NS(ManagedBuffer_get_object)( pbuffer, be_index, slot_size ) );
}

#if !defined( _GPUCODE )

SIXTRL_INLINE SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole) const*
NS(BeamElements_buffer_get_const_multipole)(
SIXTRL_BUFFER_ARGPTR_DEC const NS(Buffer) *const SIXTRL_RESTRICT buffer,
NS(buffer_size_t) const be_index )
{
typedef SIXTRL_BUFFER_DATAPTR_DEC unsigned char const* ptr_raw_t;
return NS(BeamElements_managed_buffer_get_const_multipole)(
( ptr_raw_t )( uintptr_t )NS(Buffer_get_data_begin_addr)( buffer ),
be_index, NS(Buffer_get_slot_size)( buffer ) );
}

SIXTRL_INLINE SIXTRL_BUFFER_DATAPTR_DEC NS(MultiPole)*
NS(BeamElements_buffer_get_multipole)(
SIXTRL_BUFFER_ARGPTR_DEC NS(Buffer)* SIXTRL_RESTRICT buffer,
NS(buffer_size_t) const be_index )
{
typedef SIXTRL_BUFFER_DATAPTR_DEC unsigned char* ptr_raw_t;
return NS(BeamElements_managed_buffer_get_multipole)(
( ptr_raw_t )( uintptr_t )NS(Buffer_get_data_begin_addr)( buffer ),
be_index, NS(Buffer_get_slot_size)( buffer ) );
}

SIXTRL_INLINE bool NS(MultiPole_can_be_added)(
SIXTRL_BE_ARGPTR_DEC const NS(Buffer) *const SIXTRL_RESTRICT buffer,
NS(multipole_order_t) const order,
Expand Down
16 changes: 16 additions & 0 deletions sixtracklib/common/be_multipole/be_multipole.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -364,6 +364,22 @@ namespace SIXTRL_CXX_NAMESPACE
::NS(multipole_real_t) const length = ::NS(multipole_real_t){ 0.0 },
::NS(multipole_real_t) const hxl = ::NS(multipole_real_t){ 0.0 },
::NS(multipole_real_t) const hyl = ::NS(multipole_real_t){ 0.0 } );

template<> struct ObjectTypeTraits< ::NS(MultiPole) >
{
SIXTRL_STATIC SIXTRL_INLINE object_type_id_t Type() SIXTRL_NOEXCEPT
{
return NS(OBJECT_TYPE_MULTIPOLE);
}
};

template<> struct ObjectTypeTraits< MultiPole >
{
SIXTRL_STATIC SIXTRL_INLINE object_type_id_t Type() SIXTRL_NOEXCEPT
{
return NS(OBJECT_TYPE_MULTIPOLE);
}
};
}

/* ************************************************************************* *
Expand Down
3 changes: 2 additions & 1 deletion sixtracklib/common/track/track.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -624,7 +624,8 @@ SIXTRL_INLINE NS(track_status_t) NS(Track_particle_line_objs)(
bool const finish_turn )
{
NS(track_status_t) success = SIXTRL_TRACK_SUCCESS;
bool continue_tracking = ( line_it != line_end );
bool continue_tracking = ( ( line_it != line_end ) &&
( NS(Particles_is_not_lost_value)( particles, index ) ) );

SIXTRL_ASSERT( line_it != SIXTRL_NULLPTR );
SIXTRL_ASSERT( particles != SIXTRL_NULLPTR );
Expand Down
155 changes: 154 additions & 1 deletion tests/sixtracklib/common/track/test_track_job_track_line_c99.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "sixtracklib/common/track/track_job_cpu.h"

#include <algorithm>
#include <iomanip>
#include <cstddef>
#include <cstdint>
Expand All @@ -15,13 +16,17 @@
#include "sixtracklib/testlib.h"

#include "sixtracklib/common/definitions.h"
#include "sixtracklib/common/be_drift/be_drift.h"
#include "sixtracklib/common/be_multipole/be_multipole.h"
#include "sixtracklib/common/be_limit/be_limit_rect.h"
#include "sixtracklib/common/control/definitions.h"
#include "sixtracklib/common/track/definitions.h"
#include "sixtracklib/common/buffer.h"
#include "sixtracklib/common/particles.h"
#include "sixtracklib/common/be_monitor/be_monitor.h"
#include "sixtracklib/common/track/track.h"

TEST( C99_Cpu_CpuTrackJobTrackLineTests, TrackUntilSingleParticleSetSimpleTest )
TEST( C99CpuCpuTrackJobTrackLineTests, CmpWithTrackUntilTest )
{
namespace st = SIXTRL_CXX_NAMESPACE;

Expand Down Expand Up @@ -228,4 +233,152 @@ TEST( C99_Cpu_CpuTrackJobTrackLineTests, TrackUntilSingleParticleSetSimpleTest )
::NS(Buffer_delete)( in_particles );
}


TEST( C99CpuCpuTrackJobTrackLineTests, LostParticleBehaviour )
{
using npart_t = ::NS(particle_num_elements_t);

::NS(Buffer)* elements = ::NS(Buffer_new)( ::NS(buffer_size_t){ 0 } );

::NS(DriftExact)* drift0 = ::NS(DriftExact_new)( elements );
::NS(MultiPole)* quad = ::NS(MultiPole_new)( elements, 2 );
::NS(LimitRect)* limit = ::NS(LimitRect_new)( elements );
::NS(DriftExact)* drift1 = ::NS(DriftExact_new)( elements );

drift0 = ::NS(BeamElements_buffer_get_drift_exact)( elements, 0 );
quad = ::NS(BeamElements_buffer_get_multipole)( elements, 1 );
limit = ::NS(BeamElements_buffer_get_limit_rect)( elements, 2 );
drift1 = ::NS(BeamElements_buffer_get_drift_exact)( elements, 3 );

SIXTRL_ASSERT( drift0 != nullptr );
SIXTRL_ASSERT( quad != nullptr );
SIXTRL_ASSERT( limit != nullptr );
SIXTRL_ASSERT( drift1 != nullptr );

::NS(DriftExact_set_length)( drift0, double{ 5 } );

::NS(MultiPole_set_knl_value)( quad, 0, double{ 0 } );
::NS(MultiPole_set_knl_value)( quad, 1, double{ 1e-2 } );
::NS(MultiPole_set_length)( quad, double{ 1 } );

::NS(LimitRect_set_min_x)( limit, double{ -0.1 } );
::NS(LimitRect_set_max_x)( limit, double{ 0.1 } );
::NS(LimitRect_set_min_y)( limit, double{ -0.1 } );
::NS(LimitRect_set_max_y)( limit, double{ 0.1 } );

::NS(DriftExact_set_length)( drift1, double{ 5 } );

::NS(Buffer)* cmp_pbuffer = ::NS(Buffer_new)( ::NS(buffer_size_t){ 0 } );
::NS(Buffer)* track_pbuffer = ::NS(Buffer_new)( ::NS(buffer_size_t){ 0 } );

constexpr npart_t NUM_PARTICLES = npart_t{ 100 };

::NS(Particles)* cmp_particles =
::NS(Particles_new)( cmp_pbuffer, NUM_PARTICLES );

::NS(Particles)* track_particles =
::NS(Particles_new)( track_pbuffer, NUM_PARTICLES );

double const pc0 = double{ 10e9 };
double const mass0 = double{ 1e9 };
double const q0 = double{ 1 };

double const energy0 = std::sqrt( mass0 * mass0 + pc0 * pc0 );
double const gamma0 = energy0 / mass0;
double const beta0 = std::sqrt( double{ 1 } -
double{ 1 } / ( gamma0 * gamma0 ) );

for( npart_t ii = npart_t{ 0 } ; ii < NUM_PARTICLES ; ++ii )
{
::NS(Particles_set_q0_value)( cmp_particles, ii, q0 );
::NS(Particles_set_mass0_value)( cmp_particles, ii, mass0 );
::NS(Particles_set_beta0_value)( cmp_particles, ii, beta0 );
::NS(Particles_set_gamma0_value)( cmp_particles, ii, gamma0 );
::NS(Particles_set_p0c_value)( cmp_particles, ii, pc0 );
::NS(Particles_set_x_value)(
cmp_particles, ii, static_cast< double >( ii * 0.002 ) );

::NS(Particles_set_rvv_value)( cmp_particles, ii, 1.0 );
::NS(Particles_set_rpp_value)( cmp_particles, ii, 1.0);
::NS(Particles_set_chi_value)( cmp_particles, ii, 1.0 );
::NS(Particles_set_charge_ratio_value)( cmp_particles, ii, 1.0 );

::NS(Particles_set_particle_id_value)( cmp_particles, ii, ii );
::NS(Particles_set_state_value)( cmp_particles, ii, 1 );
}

::NS(Particles_copy)( track_particles, cmp_particles );

::NS(track_status_t) track_status = ::NS(Track_all_particles_until_turn)(
cmp_particles, elements, 2 );

SIXTRL_ASSERT( track_status == ::NS(TRACK_SUCCESS) );


/* ********************************************************************* */
/* Start track_line test */

::NS(CpuTrackJob)* job = ::NS(CpuTrackJob_new)( track_pbuffer, elements );
track_status = ::NS(TrackJobNew_track_line)( job, 0, 1, false );
ASSERT_TRUE( track_status == ::NS(TRACK_SUCCESS) );

auto states_begin = ::NS(Particles_get_const_state)( track_particles );
auto states_end = states_begin;
std::advance( states_end, NUM_PARTICLES );

::NS(TrackJobNew_collect_particles)( job );

/* Check all states are equal to 1 after the first drift */
ASSERT_TRUE( std::find_if( states_begin, states_end,
[]( ::NS(particle_index_t) x ){ return x != 1; } ) == states_end );

::NS(buffer_size_t) const num_beam_elements =
::NS(Buffer_get_num_of_objects)( elements );

/* Track until the end of turn */
track_status = ::NS(TrackJobNew_track_line)(
job, 1u, num_beam_elements, true );

ASSERT_TRUE( track_status == ::NS(TRACK_SUCCESS) );
::NS(TrackJobNew_collect_particles)( job );

/* Verify that we now have lost particles due to the aperture check */
ASSERT_TRUE( std::find_if( states_begin, states_end,
[]( ::NS(particle_index_t) x ){ return x != 1; } ) != states_end );

std::vector< ::NS(particle_index_t) > const saved_states_after_first_turn(
states_begin, states_end );

SIXTRL_ASSERT( std::equal( states_begin, states_end,
saved_states_after_first_turn.begin() ) );

/* track over the first drift for the second turn */
track_status = ::NS(TrackJobNew_track_line)( job, 0u, 1u, false );
ASSERT_TRUE( track_status == ::NS(TRACK_SUCCESS) );
::NS(TrackJobNew_collect_particles)( job );

/* Since no apertures have been encountered and the global aperture
* limit should be much much larger than the current x/y values,
* the states should not have changed compared to the end of turn 1 */

ASSERT_TRUE( std::equal( states_begin, states_end,
saved_states_after_first_turn.begin() ) );

/* Finish second turn */
track_status = ::NS(TrackJobNew_track_line)(
job, 1u, num_beam_elements, true );

ASSERT_TRUE( track_status == ::NS(TRACK_SUCCESS) );
::NS(TrackJobNew_collect_particles)( job );

/* Compare against the results obtained by performing
* NS(Track_all_particles_until_turn) for the whole two turns
* in one go */

double const ABS_TOLERANCE = double{ 1e-14 };

ASSERT_TRUE( ::NS(Particles_compare_real_values_with_treshold)(
cmp_particles, track_particles, ABS_TOLERANCE ) == 0 );
}

/* end: tests/sixtracklib/common/track/test_track_job_track_until_c99.cpp */
Loading

0 comments on commit 947643d

Please sign in to comment.