forked from celeritas-project/celeritas
/
OrangeTrackView.hh
1117 lines (955 loc) · 34.9 KB
/
OrangeTrackView.hh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
//----------------------------------*-C++-*----------------------------------//
// Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file orange/OrangeTrackView.hh
//---------------------------------------------------------------------------//
#pragma once
#include "corecel/Macros.hh"
#include "corecel/Types.hh"
#include "corecel/cont/Array.hh"
#include "corecel/sys/ThreadId.hh"
#include "OrangeData.hh"
#include "OrangeTypes.hh"
#include "transform/TransformVisitor.hh"
#include "univ/SimpleUnitTracker.hh"
#include "univ/TrackerVisitor.hh"
#include "univ/UniverseTypeTraits.hh"
#include "univ/detail/Types.hh"
#include "detail/LevelStateAccessor.hh"
#include "detail/UniverseIndexer.hh"
namespace celeritas
{
//---------------------------------------------------------------------------//
/*!
* Navigate through an ORANGE geometry on a single thread.
*
* Ordering requirements:
* - initialize (through assignment) must come first
* - access (pos, dir, volume/surface/is_outside/is_on_boundary) good at any
* time
* - \c find_next_step
* - \c find_safety or \c move_internal or \c move_to_boundary
* - if on boundary, \c cross_boundary
* - at any time, \c set_dir , but then must do \c find_next_step before any
* \c move or \c cross action above
*
* The main point is that \c find_next_step depends on the current
* straight-line direction, \c move_to_boundary and \c move_internal (with
* a step length) depends on that distance, and
* \c cross_boundary depends on being on the boundary with a knowledge of the
* post-boundary state.
*
* \c move_internal with a position \em should depend on the safety distance
* but that's not yet implemented.
*/
class OrangeTrackView
{
public:
//!@{
//! \name Type aliases
using ParamsRef = NativeCRef<OrangeParamsData>;
using StateRef = NativeRef<OrangeStateData>;
using Initializer_t = GeoTrackInitializer;
//!@}
//! Helper struct for initializing from an existing geometry state
struct DetailedInitializer
{
OrangeTrackView const& other; //!< Existing geometry
Real3 const& dir; //!< New direction
};
public:
// Construct from params and state
inline CELER_FUNCTION OrangeTrackView(ParamsRef const& params,
StateRef const& states,
TrackSlotId tid);
// Initialize the state
inline CELER_FUNCTION OrangeTrackView& operator=(Initializer_t const& init);
// Initialize the state from a parent state and new direction
inline CELER_FUNCTION OrangeTrackView&
operator=(DetailedInitializer const& init);
//// ACCESSORS ////
// The current position
inline CELER_FUNCTION Real3 const& pos() const;
// The current direction
inline CELER_FUNCTION Real3 const& dir() const;
// The current volume ID (null if outside)
inline CELER_FUNCTION VolumeId volume_id() const;
// The current surface ID
inline CELER_FUNCTION SurfaceId surface_id() const;
// After 'find_next_step', the next straight-line surface
inline CELER_FUNCTION SurfaceId next_surface_id() const;
// Whether the track is outside the valid geometry region
inline CELER_FUNCTION bool is_outside() const;
// Whether the track is exactly on a surface
inline CELER_FUNCTION bool is_on_boundary() const;
//// OPERATIONS ////
// Find the distance to the next boundary
inline CELER_FUNCTION Propagation find_next_step();
// Find the distance to the next boundary, up to and including a step
inline CELER_FUNCTION Propagation find_next_step(real_type max_step);
// Find the distance to the nearest boundary in any direction
inline CELER_FUNCTION real_type find_safety();
// Find the distance to the nearest nearby boundary in any direction
inline CELER_FUNCTION real_type find_safety(real_type max_step);
// Move to the boundary in preparation for crossing it
inline CELER_FUNCTION void move_to_boundary();
// Move within the volume
inline CELER_FUNCTION void move_internal(real_type step);
// Move within the volume to a specific point
inline CELER_FUNCTION void move_internal(Real3 const& pos);
// Cross from one side of the current surface to the other
inline CELER_FUNCTION void cross_boundary();
// Change direction
inline CELER_FUNCTION void set_dir(Real3 const& newdir);
private:
//// TYPES ////
using LSA = detail::LevelStateAccessor;
//// DATA ////
ParamsRef const& params_;
StateRef const& states_;
TrackSlotId track_slot_;
//// PRIVATE STATE MUTATORS ////
// The current level
inline CELER_FUNCTION void level(LevelId);
// The boundary on the current surface level
inline CELER_FUNCTION void boundary(BoundaryResult);
// The next step distance, as stored on the state
inline CELER_FUNCTION void next_step(real_type dist);
// The next surface to be encounted
inline CELER_FUNCTION void next_surf(detail::OnLocalSurface const&);
// The level of the next surface to be encounted
inline CELER_FUNCTION void next_surface_level(LevelId);
//// PRIVATE STATE ACCESSORS ////
// The current level
inline CELER_FUNCTION LevelId const& level() const;
// The current surface level
inline CELER_FUNCTION LevelId const& surface_level() const;
// The local surface on the current surface level
inline CELER_FUNCTION LocalSurfaceId const& surf() const;
// The sense on the current surface level
inline CELER_FUNCTION Sense const& sense() const;
// The boundary on the current surface level
inline CELER_FUNCTION BoundaryResult const& boundary() const;
// The next step distance, as stored on the state
inline CELER_FUNCTION real_type const& next_step() const;
// The next surface to be encounted
inline CELER_FUNCTION detail::OnLocalSurface next_surf() const;
// The level of the next surface to be encounted
inline CELER_FUNCTION LevelId const& next_surface_level() const;
//// HELPER FUNCTIONS ////
// Iterate over lower levels to find the next step
inline CELER_FUNCTION Propagation
find_next_step_impl(detail::Intersection isect);
// Create local sense reference
inline CELER_FUNCTION Span<Sense> make_temp_sense() const;
// Create local distance
inline CELER_FUNCTION detail::TempNextFace make_temp_next() const;
inline CELER_FUNCTION detail::LocalState
make_local_state(LevelId level) const;
// Whether the next distance-to-boundary has been found
inline CELER_FUNCTION bool has_next_step() const;
// Whether the next surface has been found
inline CELER_FUNCTION bool has_next_surface() const;
// Invalidate the next distance-to-boundary and surface
inline CELER_FUNCTION void clear_next();
// Assign the surface on the current level
inline CELER_FUNCTION void
surface(LevelId level, detail::OnLocalSurface surf);
// Clear the surface on the current level
inline CELER_FUNCTION void clear_surface();
// Make a LevelStateAccessor for the current thread and level
inline CELER_FUNCTION LSA make_lsa() const;
// Make a LevelStateAccessor for the current thread and a given level
inline CELER_FUNCTION LSA make_lsa(LevelId level) const;
// Get the daughter ID for the volume in the universe (or null)
inline CELER_FUNCTION DaughterId get_daughter(LSA const& lsa);
// Get the transform ID for the given daughter.
inline CELER_FUNCTION TransformId get_transform(DaughterId daughter_id);
// Get the transform ID to move from this level to the one below.
inline CELER_FUNCTION TransformId get_transform(LevelId lev);
};
//---------------------------------------------------------------------------//
// MEMBER FUNCTIONS
//---------------------------------------------------------------------------//
/*!
* Construct from persistent and state data.
*/
CELER_FUNCTION
OrangeTrackView::OrangeTrackView(ParamsRef const& params,
StateRef const& states,
TrackSlotId tid)
: params_(params), states_(states), track_slot_(tid)
{
CELER_EXPECT(params_);
CELER_EXPECT(states_);
CELER_EXPECT(track_slot_ < states.size());
}
//---------------------------------------------------------------------------//
/*!
* Construct the state.
*
* Expensive. This function should only be called to initialize an event from a
* starting location and direction. Secondaries will initialize their states
* from a copy of the parent.
*/
CELER_FUNCTION OrangeTrackView&
OrangeTrackView::operator=(Initializer_t const& init)
{
CELER_EXPECT(is_soft_unit_vector(init.dir));
// Create local state
detail::LocalState local;
local.pos = init.pos;
local.dir = init.dir;
local.volume = {};
local.surface = {};
local.temp_sense = this->make_temp_sense();
// Helpers for applying parent-to-daughter transformations
TransformVisitor apply_transform{params_};
auto transform_down_local = [&local](auto&& t) {
local.pos = t.transform_down(local.pos);
local.dir = t.rotate_down(local.dir);
};
// Recurse into daughter universes starting with the outermost universe
UniverseId uid = top_universe_id();
DaughterId daughter_id;
size_type level = 0;
do
{
TrackerVisitor visit_tracker{params_};
auto tinit = visit_tracker(
[&local](auto&& t) { return t.initialize(local); }, uid);
// TODO: error correction/graceful failure if initialization failed
CELER_ASSERT(tinit.volume && !tinit.surface);
auto lsa = this->make_lsa(LevelId{level});
lsa.vol() = tinit.volume;
lsa.pos() = local.pos;
lsa.dir() = local.dir;
lsa.universe() = uid;
CELER_ASSERT(tinit.volume);
daughter_id = visit_tracker(
[&tinit](auto&& t) { return t.daughter(tinit.volume); }, uid);
if (daughter_id)
{
auto const& daughter = params_.daughters[daughter_id];
// Apply "transform down" based on stored transform
apply_transform(transform_down_local, daughter.transform_id);
// Update universe and increase level depth
uid = daughter.universe_id;
++level;
}
} while (daughter_id);
// Save fiound level
this->level(LevelId{level});
// Set default boundary condition
this->boundary(BoundaryResult::exiting);
// Clear surface information
this->clear_surface();
this->clear_next();
CELER_ENSURE(!this->has_next_step());
return *this;
}
//---------------------------------------------------------------------------//
/*!
* Construct the state from a direction and a copy of the parent state.
*/
CELER_FUNCTION
OrangeTrackView& OrangeTrackView::operator=(DetailedInitializer const& init)
{
CELER_EXPECT(is_soft_unit_vector(init.dir));
if (this != &init.other)
{
// Copy init track's position and logical state
this->level(states_.level[init.other.track_slot_]);
this->surface(init.other.surface_level(),
{init.other.surf(), init.other.sense()});
this->boundary(init.other.boundary());
for (auto lev : range(LevelId{this->level() + 1}))
{
// Copy all data accessed via LSA
auto lsa = this->make_lsa(lev);
lsa = init.other.make_lsa(lev);
}
}
// Clear the next step information since we're changing direction or
// initializing a new state
this->clear_next();
// Transform direction from global to local
Real3 localdir = init.dir;
auto apply_transform = TransformVisitor{params_};
auto rotate_down
= [&localdir](auto&& t) { localdir = t.rotate_down(localdir); };
for (auto lev : range(LevelId{this->level()}))
{
auto lsa = this->make_lsa(lev);
lsa.dir() = localdir;
apply_transform(rotate_down,
this->get_transform(this->get_daughter(lsa)));
}
// Save final level
this->make_lsa().dir() = localdir;
CELER_ENSURE(!this->has_next_step());
return *this;
}
//---------------------------------------------------------------------------//
/*!
* The current position.
*/
CELER_FUNCTION Real3 const& OrangeTrackView::pos() const
{
return this->make_lsa(LevelId{0}).pos();
}
//---------------------------------------------------------------------------//
/*!
* The current direction.
*/
CELER_FUNCTION Real3 const& OrangeTrackView::dir() const
{
return this->make_lsa(LevelId{0}).dir();
}
//---------------------------------------------------------------------------//
/*!
* The current volume ID.
*
* \note It is allowable to call this function when "outside", because the
* outside in ORANGE is just a special volume. Other geometries may not have
* that behavior.
*/
CELER_FUNCTION VolumeId OrangeTrackView::volume_id() const
{
auto lsa = this->make_lsa();
detail::UniverseIndexer ui(params_.universe_indexer_data);
return ui.global_volume(lsa.universe(), lsa.vol());
}
//---------------------------------------------------------------------------//
/*!
* The current surface ID.
*/
CELER_FUNCTION SurfaceId OrangeTrackView::surface_id() const
{
if (this->is_on_boundary())
{
auto lsa = this->make_lsa(this->surface_level());
detail::UniverseIndexer ui{params_.universe_indexer_data};
return ui.global_surface(lsa.universe(), this->surf());
}
else
{
return SurfaceId{};
}
}
//---------------------------------------------------------------------------//
/*!
* After 'find_next_step', the next straight-line surface.
*/
CELER_FUNCTION SurfaceId OrangeTrackView::next_surface_id() const
{
CELER_EXPECT(this->has_next_surface());
auto lsa = this->make_lsa(this->next_surface_level());
detail::UniverseIndexer ui{params_.universe_indexer_data};
return ui.global_surface(lsa.universe(), this->next_surf().id());
}
//---------------------------------------------------------------------------//
/*!
* Whether the track is outside the valid geometry region.
*/
CELER_FUNCTION bool OrangeTrackView::is_outside() const
{
// Zeroth volume in outermost universe is always the exterior by
// construction in ORANGE
auto lsa = this->make_lsa(LevelId{0});
return lsa.vol() == orange_exterior_volume;
}
//---------------------------------------------------------------------------//
/*!
* Whether the track is exactly on a surface.
*/
CELER_FORCEINLINE_FUNCTION bool OrangeTrackView::is_on_boundary() const
{
return static_cast<bool>(this->surface_level());
}
//---------------------------------------------------------------------------//
/*!
* Find the distance to the next geometric boundary.
*/
CELER_FUNCTION Propagation OrangeTrackView::find_next_step()
{
if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
{
// On a boundary, headed back in: next step is zero
return {0, true};
}
// Find intersection at the top level: always the first simple unit
auto global_isect = [this] {
SimpleUnitTracker t{params_, SimpleUnitId{0}};
return t.intersect(this->make_local_state(LevelId{0}));
}();
// Find intersection for all lower levels
return this->find_next_step_impl(global_isect);
}
//---------------------------------------------------------------------------//
/*!
* Find a nearby distance to the next geometric boundary up to a distance.
*
* This may reduce the number of surfaces needed to check, sort, or write to
* temporary memory, thereby speeding up transport.
*/
CELER_FUNCTION Propagation OrangeTrackView::find_next_step(real_type max_step)
{
CELER_EXPECT(max_step > 0);
if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
{
// On a boundary, headed back in: next step is zero
return {0, true};
}
// Find intersection at the top level: always the first simple unit
auto global_isect = [this, &max_step] {
SimpleUnitTracker t{params_, SimpleUnitId{0}};
return t.intersect(this->make_local_state(LevelId{0}), max_step);
}();
// Find intersection for all lower levels
auto result = this->find_next_step_impl(global_isect);
CELER_ENSURE(result.distance <= max_step);
return result;
}
//---------------------------------------------------------------------------//
/*!
* Move to the next straight-line boundary but do not change volume.
*/
CELER_FUNCTION void OrangeTrackView::move_to_boundary()
{
CELER_EXPECT(this->boundary() != BoundaryResult::reentrant);
CELER_EXPECT(this->has_next_step());
CELER_EXPECT(this->has_next_surface());
// Physically move next step
real_type const dist = this->next_step();
for (auto i : range(this->level() + 1))
{
auto lsa = this->make_lsa(LevelId{i});
axpy(dist, lsa.dir(), &lsa.pos());
}
this->surface(this->next_surface_level(), this->next_surf());
this->clear_next();
CELER_ENSURE(this->is_on_boundary());
}
//---------------------------------------------------------------------------//
/*!
* Move within the current volume.
*
* The straight-line distance *must* be less than the distance to the
* boundary.
*/
CELER_FUNCTION void OrangeTrackView::move_internal(real_type dist)
{
CELER_EXPECT(this->has_next_step());
CELER_EXPECT(dist > 0 && dist <= this->next_step());
CELER_EXPECT(dist != this->next_step() || !this->has_next_surface());
// Move and update the next step
for (auto i : range(this->level() + 1))
{
auto lsa = this->make_lsa(LevelId{i});
axpy(dist, lsa.dir(), &lsa.pos());
}
this->next_step(this->next_step() - dist);
this->clear_surface();
}
//---------------------------------------------------------------------------//
/*!
* Move within the current volume to a nearby point.
*
* \todo Currently it's up to the caller to make sure that the position is
* "nearby". We should actually test this with an "is inside" call.
*/
CELER_FUNCTION void OrangeTrackView::move_internal(Real3 const& pos)
{
// Transform all nonlocal levels
auto local_pos = pos;
auto apply_transform = TransformVisitor{params_};
auto translate_down
= [&local_pos](auto&& t) { local_pos = t.transform_down(local_pos); };
for (auto lev : range(LevelId{this->level()}))
{
auto lsa = this->make_lsa(lev);
lsa.pos() = local_pos;
// Apply "transform down" based on stored transform
apply_transform(translate_down,
this->get_transform(this->get_daughter(lsa)));
}
// Save final level
auto lsa = this->make_lsa();
lsa.pos() = local_pos;
// Clear surface state and next-step info
this->clear_surface();
this->clear_next();
}
//---------------------------------------------------------------------------//
/*!
* Cross from one side of the current surface to the other.
*
* The position *must* be on the boundary following a move-to-boundary. This
* should only be called once per boundary crossing.
*/
CELER_FUNCTION void OrangeTrackView::cross_boundary()
{
CELER_EXPECT(this->is_on_boundary());
CELER_EXPECT(!this->has_next_step());
auto sl = this->surface_level();
auto lsa = this->make_lsa(sl);
if (CELER_UNLIKELY(this->boundary() == BoundaryResult::reentrant))
{
// Direction changed while on boundary leading to no change in
// volume/surface. This is logically equivalent to a reflection.
this->boundary(BoundaryResult::exiting);
return;
}
// Flip current sense from "before crossing" to "after"
detail::LocalState local;
local.pos = lsa.pos();
local.dir = lsa.dir();
local.volume = lsa.vol();
local.surface = {this->surf(), flip_sense(this->sense())};
local.temp_sense = this->make_temp_sense();
// Helpers for updating local transforms
TransformVisitor apply_transform{params_};
auto transform_down_local = [&local](auto&& t) {
local.pos = t.transform_down(local.pos);
local.dir = t.rotate_down(local.dir);
};
// Update the post-crossing volume
TrackerVisitor visit_tracker{params_};
auto tinit = visit_tracker(
[&local](auto&& t) { return t.cross_boundary(local); }, lsa.universe());
CELER_ASSERT(tinit.volume);
if (!CELERITAS_DEBUG && CELER_UNLIKELY(!tinit.volume))
{
// Initialization failure on release mode: set to exterior volume
// rather than segfaulting
// TODO: error correction or more graceful failure than losing
// energy
tinit.volume = orange_exterior_volume;
tinit.surface = {};
}
lsa.vol() = tinit.volume;
this->surface(sl, tinit.surface);
this->boundary(BoundaryResult::exiting);
// Starting with the current level (i.e., next_surface_level), iterate
// down into the deepest level
size_type level = sl.get();
LocalVolumeId volume = tinit.volume;
auto universe_id = lsa.universe();
CELER_ASSERT(volume);
auto daughter_id = visit_tracker(
[&volume](auto&& t) { return t.daughter(volume); }, lsa.universe());
while (daughter_id)
{
auto const& daughter = params_.daughters[daughter_id];
// Make the current level the daughter level
++level;
universe_id = daughter.universe_id;
// Create local state on the daughter level
apply_transform(transform_down_local, daughter.transform_id);
local.volume = {};
local.surface = {};
// Initialize in daughter and get IDs of volume and potential daughter
volume = visit_tracker(
[&local](auto&& t) { return t.initialize(local).volume; },
universe_id);
CELER_ASSERT(volume);
daughter_id = visit_tracker(
[&volume](auto&& t) { return t.daughter(volume); }, universe_id);
auto lsa = make_lsa(LevelId{level});
lsa.vol() = volume;
lsa.pos() = local.pos;
lsa.dir() = local.dir;
lsa.universe() = universe_id;
}
this->level(LevelId{level});
CELER_ENSURE(this->is_on_boundary());
}
//---------------------------------------------------------------------------//
/*!
* Change the track's direction.
*
* This happens after a scattering event or movement inside a magnetic field.
* It resets the calculated distance-to-boundary. It is allowed to happen on
* the boundary, but changing direction so that it goes from pointing outward
* to inward (or vice versa) will mean that \c cross_boundary will be a
* null-op.
*
* TODO: This needs to be updated to handle reflections through levels
*/
CELER_FUNCTION void OrangeTrackView::set_dir(Real3 const& newdir)
{
CELER_EXPECT(is_soft_unit_vector(newdir));
if (this->is_on_boundary())
{
// Changing direction on a boundary, which may result in not leaving
// current volume upon the cross_surface call
auto lsa = this->make_lsa(this->surface_level());
TrackerVisitor visit_tracker{params_};
auto normal = visit_tracker(
[pos = lsa.pos(), local_surface = this->surf()](auto&& t) {
return t.normal(pos, local_surface);
},
lsa.universe());
// Normal is in *local* coordinates but newdir is in *global*: rotate
// up to check
auto apply_transform = TransformVisitor{params_};
auto rotate_up = [&normal](auto&& t) { normal = t.rotate_up(normal); };
for (auto level : range<int>(this->level().unchecked_get()).step(-1))
{
apply_transform(rotate_up, this->get_transform(LevelId(level)));
}
// Evaluate whether the direction dotted with the surface normal
// changes (i.e. heading from inside to outside or vice versa).
if ((dot_product(normal, newdir) >= 0)
!= (dot_product(normal, this->dir()) >= 0))
{
// The boundary crossing direction has changed! Reverse our
// plans to change the logical state and move to a new volume.
this->boundary(flip_boundary(this->boundary()));
}
}
// Complete direction setting by transforming direction all the way down
Real3 localdir = newdir;
auto apply_transform = TransformVisitor{params_};
auto rotate_down
= [&localdir](auto&& t) { localdir = t.rotate_down(localdir); };
for (auto lev : range(this->level()))
{
auto lsa = this->make_lsa(lev);
lsa.dir() = localdir;
apply_transform(rotate_down,
this->get_transform(this->get_daughter(lsa)));
}
// Save final level
this->make_lsa().dir() = localdir;
this->clear_next();
}
//---------------------------------------------------------------------------//
// STATE ACCESSORS
//---------------------------------------------------------------------------//
/*!
* The current level.
*/
CELER_FORCEINLINE_FUNCTION void OrangeTrackView::level(LevelId lev)
{
states_.level[track_slot_] = lev;
}
/*!
* The boundary on the current surface level.
*/
CELER_FORCEINLINE_FUNCTION void OrangeTrackView::boundary(BoundaryResult br)
{
states_.boundary[track_slot_] = br;
}
/*!
* The next step distance.
*/
CELER_FORCEINLINE_FUNCTION void OrangeTrackView::next_step(real_type dist)
{
states_.next_step[track_slot_] = dist;
}
/*!
* The next surface to be encountered.
*/
CELER_FORCEINLINE_FUNCTION void
OrangeTrackView::next_surf(detail::OnLocalSurface const& s)
{
states_.next_surf[track_slot_] = s.id();
states_.next_sense[track_slot_] = s.unchecked_sense();
}
/*!
* The level of the next surface to be encounted.
*/
CELER_FORCEINLINE_FUNCTION void OrangeTrackView::next_surface_level(LevelId lev)
{
states_.next_level[track_slot_] = lev;
}
//---------------------------------------------------------------------------//
// CONST STATE ACCESSORS
//---------------------------------------------------------------------------//
/*!
* The current level.
*/
CELER_FORCEINLINE_FUNCTION LevelId const& OrangeTrackView::level() const
{
return states_.level[track_slot_];
}
/*!
* The current surface level.
*/
CELER_FORCEINLINE_FUNCTION LevelId const& OrangeTrackView::surface_level() const
{
return states_.surface_level[track_slot_];
}
/*!
* The local surface on the current surface level.
*/
CELER_FORCEINLINE_FUNCTION LocalSurfaceId const& OrangeTrackView::surf() const
{
return states_.surf[track_slot_];
}
/*!
* The sense on the current surface level.
*/
CELER_FORCEINLINE_FUNCTION Sense const& OrangeTrackView::sense() const
{
return states_.sense[track_slot_];
}
/*!
* The boundary on the current surface level.
*/
CELER_FORCEINLINE_FUNCTION BoundaryResult const&
OrangeTrackView::boundary() const
{
return states_.boundary[track_slot_];
}
/*!
* The next step distance.
*/
CELER_FORCEINLINE_FUNCTION real_type const& OrangeTrackView::next_step() const
{
return states_.next_step[track_slot_];
}
/*!
* The next surface to be encountered.
*/
CELER_FORCEINLINE_FUNCTION detail::OnLocalSurface
OrangeTrackView::next_surf() const
{
return {states_.next_surf[track_slot_], states_.next_sense[track_slot_]};
}
/*!
* The level of the next surface to be encounted.
*/
CELER_FORCEINLINE_FUNCTION LevelId const&
OrangeTrackView::next_surface_level() const
{
return states_.next_level[track_slot_];
}
//---------------------------------------------------------------------------//
// HELPER FUNCTIONS
//---------------------------------------------------------------------------//
/*!
* Iterate over levels 1 to N to find the next step.
*
* Caller is responsible for finding the candidate next step on level 0, and
* passing the resultant Intersection object as an argument.
*/
CELER_FUNCTION Propagation
OrangeTrackView::find_next_step_impl(detail::Intersection isect)
{
TrackerVisitor visit_tracker{params_};
// The level with minimum distance to intersection
LevelId min_level{0};
// Find the nearest intersection from level 0 to current level
// inclusive, prefering the shallowest level (i.e., lowest uid)
for (auto levelid : range(LevelId{1}, this->level() + 1))
{
auto uid = this->make_lsa(levelid).universe();
auto local_isect = visit_tracker(
[local_state = this->make_local_state(levelid), &isect](auto&& t) {
return t.intersect(local_state, isect.distance);
},
uid);
if (local_isect.distance < isect.distance)
{
isect = local_isect;
min_level = levelid;
}
}
this->next_step(isect.distance);
this->next_surf(isect.surface);
if (isect)
{
// Save level corresponding to the intersection
this->next_surface_level(min_level);
}
Propagation result;
result.distance = isect.distance;
result.boundary = static_cast<bool>(isect);
return result;
}
//---------------------------------------------------------------------------//
/*!
* Find the distance to the nearest boundary in any direction.
*
* The safety distance at a given point is the minimum safety distance over all
* levels, since surface deduplication can potentionally elide bounding
* surfaces at more deeply embedded levels.
*/
CELER_FUNCTION real_type OrangeTrackView::find_safety()
{
CELER_EXPECT(!this->is_on_boundary());
TrackerVisitor visit_tracker{params_};
real_type min_safety_dist = numeric_limits<real_type>::infinity();
for (auto lev : range(LevelId{this->level() + 1}))
{
auto lsa = this->make_lsa(lev);
auto sd = visit_tracker(
[&lsa](auto&& t) { return t.safety(lsa.pos(), lsa.vol()); },
lsa.universe());
min_safety_dist = celeritas::min(min_safety_dist, sd);
}
return min_safety_dist;
}
//---------------------------------------------------------------------------//
/*!
* Find the distance to the nearest nearby boundary.
*
* Since we currently support only "simple" safety distances, we can't
* eliminate anything by checking only nearby surfaces.
*/
CELER_FUNCTION real_type OrangeTrackView::find_safety(real_type)
{
return this->find_safety();
}
//---------------------------------------------------------------------------//
/*!
* Get a reference to the current volume, or to world volume if outside.
*/
CELER_FUNCTION Span<Sense> OrangeTrackView::make_temp_sense() const
{
auto const max_faces = params_.scalars.max_faces;
auto offset = track_slot_.get() * max_faces;
return states_.temp_sense[AllItems<Sense, MemSpace::native>{}].subspan(
offset, max_faces);
}
//---------------------------------------------------------------------------//
/*!
* Set up intersection scratch space.
*/
CELER_FUNCTION detail::TempNextFace OrangeTrackView::make_temp_next() const
{
auto const max_isect = params_.scalars.max_intersections;
auto offset = track_slot_.get() * max_isect;
detail::TempNextFace result;
result.face = states_.temp_face[AllItems<FaceId>{}].data() + offset;
result.distance = states_.temp_distance[AllItems<real_type>{}].data()
+ offset;
result.isect = states_.temp_isect[AllItems<size_type>{}].data() + offset;
result.size = max_isect;
return result;
}
//---------------------------------------------------------------------------//
/*!
* Create a local state.
*/
CELER_FUNCTION detail::LocalState
OrangeTrackView::make_local_state(LevelId level) const
{
detail::LocalState local;
auto lsa = this->make_lsa(level);
local.pos = lsa.pos();
local.dir = lsa.dir();
local.volume = lsa.vol();
if (level == this->surface_level())
{