* Bialek et al. [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak.  Statistical mechanics for natural flocks of birds.  _Proceedings of the National Academy of Sciences_ , 109(13):4786–4791, May 2012.  ISSN 1091-6490.  doi: 10.1073/pnas.1118633109.  URL http://dx.doi.org/10.1073/pnas.1118633109.
  * Chang et al. [2015] Angel X. Chang, Thomas Funkhouser, Leonidas Guibas, Pat Hanrahan, Qixing Huang, Zimo Li, Silvio Savarese, Manolis Savva, Shuran Song, Hao Su, Jianxiong Xiao, Li Yi, and Fisher Yu.  ShapeNet: An Information-Rich 3D Model Repository.  Technical Report arXiv:1512.03012 [cs.GR], Stanford University — Princeton University — Toyota Technological Institute at Chicago, 2015.
  * Menze et al. [2018] Moritz Menze, Christian Heipke, and Andreas Geiger.  Object scene flow.  _ISPRS Journal of Photogrammetry and Remote Sensing (JPRS)_ , 2018\.
  * Vassileios Balntas and Mikolajczyk [2016] Daniel Ponsa Vassileios Balntas, Edgar Riba and Krystian Mikolajczyk.  Learning local feature descriptors with triplets and shallow convolutional neural networks.  In Edwin R. Hancock Richard C. Wilson and William A. P. Smith, editors, _Proceedings of the British Machine Vision Conference (BMVC)_ , pages 119.1–119.11. BMVA Press, Sep 2016.  ISBN 1-901725-59-6.  doi: 10.5244/C.30.119.  URL https://dx.doi.org/10.5244/C.30.119.
  * Hoffer and Ailon [2015] Elad Hoffer and Nir Ailon.  Deep metric learning using triplet network.  In Yoshua Bengio and Yann LeCun, editors, _ICLR (Workshop)_ , 2015\.  URL http://dblp.uni-trier.de/db/conf/iclr/iclr2015w.html#HofferA14.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, 2016.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, editors, _Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015_ , pages 234–241, Cham, 2015. Springer International Publishing.  ISBN 978-3-319-24574-4.
  * Karnewar and Wang [2020] Animesh Karnewar and Oliver Wang.  Msg-gan: Multi-scale gradients for generative adversarial networks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2020.
  * Geiger et al. [2012] Andreas Geiger, Philip Lenz, and Raquel Urtasun.  Are we ready for autonomous driving? the kitti vision benchmark suite.  In _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2012.
  * Geiger et al. [2013] Andreas Geiger, Philip Lenz, Christoph Stiller, and Raquel Urtasun.  Vision meets robotics: The kitti dataset.  _International Journal of Robotics Research (IJRR)_ , 2013.
  * Liu et al. [2018] Xingyu Liu, Charles Ruizhongtai Qi, and Leonidas J. Guibas.  Flownet3d: Learning scene flow in 3d point clouds.  https://github.com/xingyul/flownet3d, 2018.
  * Kuhn and Yaw [1955] H. W. Kuhn and Bryn Yaw.  The hungarian method for the assignment problem.  _Naval Res. Logist. Quart_ , pages 83–97, 1955.

###### List of Figures

  1. 1.1 Illustration of the scene flow task. Blue points belong to a snapshot at time $t$ and purple points belong to a snapshot at time $t+1$. It is possible that a different number of points is recorded at different time steps and that there is no one-to-one correspondence between points of different frames. The task of scene flow is to estimate motion vectors for each point of a scene.
  2. 3.1 Illustrative example of the Correspondence and the Re-sampling Mechanisms. The blue points belong to $C_{1}$ and the purple points belong to $C_{2}$. When correspondence is assumed, each point at $C_{2}$ has a corresponding point at $C_{1}$. When Re-Sampling is assumed this correspondence is not present any longer, and the number of points may be different at each time step.
  3. 3.2 Snapshot of a flock (A) and the velocity vector of individual birds (B). Figure taken from [38].
  4. 3.3 Illustrative example of the Correspondence and the Re-sampling Mechanism. The shaded colors represent the areas of the objects that are not visible when using one or other assumption. F.o.V stands for Field of View.
  5. 4.1 Illustration of the setup for Adversarial Metric Training
  6. 4.2 Examples of sub-sampled point clouds visualized by the black and red dots. From left to right the datasets are Single ShapeNet, Multi ShapeNet, FlyingThings3D, KITTI. The points of different colors have uniform coverage of the scene and are expected to represent each object equally well.
  7. 4.3 Illustration of the inputs for the loss calculation for the Flow Extractor and for the Cloud Embedder.
  8. (a) Flow Extractor is trained to minimize the L2 distance between $C_{t+1}$ and $\hat{C}_{t+1}$ in the latent space. Flow Extractor is trained to minimize the L2 distance between $C_{t+1}$ and $\hat{C}_{t+1}$ in the latent space. Flow Extractor is trained to minimize the L2 space.
  9. (b) Cloud Embedder is trained to minimize the triplet margin between $pc_{t+1}$ and $\hat{pc}_{t+1}$. Two non-overlapping sub-sets of $pc_{t+1}$ are randomly sampled. The sampling is illustrated as masked areas of the point clouds, however the actual sampling is agnostic to any aspects of the point cloud. That means that all regions are, in expectation, equally represented in both sub-sets.
  10. 4.4 Illustration of the Multi-Scale Triplet Loss. The Feature Extractors can be any differentiable model, in our work we use PointNet++ [26] Feature Extractors.
  11. 5.1 Examples of frames from Single ShapeNet. Colors indicate the sections of the objects.
  12. (a) Scene with table.
  13. (b) Scene with airplane.
  14. 5.2 Examples of frames from Multi ShapeNet.
  15. (a) Scene with airplane (pink) and table (gray) and vase (green).
  16. (b) Scene with airplane (light purple) and table (gray).
  17. 5.3 Examples of frames from FlowNet3D. The points of 5.3(b) and 5.3(d) were colored with the RGB from 5.3(a) and 5.3(c).
  18. (a) Image of scene 5.3(b)
  19. (b) Point cloud of scene 5.3(a)
  20. (c) Image of scene 5.3(d)
  21. (d) Point cloud of scene 5.3(c)
  22. 5.4 Examples of frames from KITTI. Height is used to color the points.
  23. (a) Example of KITTI with ground.
  24. (b) Example of KITTI with ground
  25. (c) Example of KITTI without ground
  26. (d) Example of KITTI without ground
  27. 5.5 Examples of frames from Lyft. Height is used to color the points.
  28. (a)
  29. (b)
  30. 6.1 Example taken from Single ShapeNet. The 2D projection shows the Ground Truth flow vectors and the estimations made by the different baselines. The blue points correspond to time $t$ and orange points to time $t+1$.
  31. 6.2 Example taken from Multi ShapeNet. The 2D projection shows the Ground Truth flow vectors and the estimations made by the different baselines. The blue points correspond to time $t$ and orange points to time $t+1$.
  32. 6.3 Example taken from FlyingThings3D. The 2D projection is chosen to make occlusions evident. The pink dotted line highlights an object that was present in time $t=0$ and was occluded at time $t=1$. The blue points correspond to time $t$ and orange points to time $t+1$.
  33. 6.4 Examples of Segmenter and FlowNet3D on the synthetic datasets.
  34. (a) Qualitative results from Segmenter.
  35. (b) Qualitative results from FlowNet3D.
  36. 6.5 Qualitative results from FlowNet3D after domain adaptation.
  37. (a) width=0.9
  38. (b) width=0.9
  39. 6.6 Qualitative results from Adversarial Metric Learning on Single ShapeNet, Multi ShapeNet and FlyingThings3D. The pink doted lines shows examples of occlusions that are not correctly handled by the flow extractors.
  40. (a) FlowNet3D used as flow extractor. The pink doted line shows an example of occlusion.
  41. (b) PointPWC-net used as flow extractor.
  42. 6.7 Qualitative results from Adversarial Metric Learning on KITTI with and without ground.
  43. (a) FlowNet3D used as flow extractor.
  44. (b) PointPWC-net used as flow extractor.
  45. 6.8 Qualitative results from FlowNet3D on Single ShapeNet.
  46. (a) 2D view of a table from Single ShapeNet. Flow Extractor used is FlowNet3D. Blue Points belong to frame 1, orange points belong to frame 2 and the black arrows are the flow vectors.
  47. (b) 2D projection of an Airplane from Single ShapeNet. Flow Extractor cluster points in order to minimize the nearest neighbor distance. The top left image shows the target point cloud, the top right image shows the predicted one. The flow collapses most points around the centroid of the object. On bottom left the ground truth flow vectors are shown, on the bottom right the estimated flow vectors are displayed.
  48. 6.9 Illustrative example of the different set distances. Blue point belong to frame $t$, purple points belong to frame $t+1$. Gray arrows represent the ground truth flow. The blue arrows represent the pointwise distances with respect to $t$ and the purple the point wise distances with respect to $t+1$. The set distances are the sum of the arrows in each image. The image makes explicit how the set distances are ill-suited when the Re-sampling mechanism is used. Nearest neighbor tends to underestimate the total motion in both mechanisms. Chamfer distance has the same underestimation problem, however the bilateral assignment regularizes it to some extent. Hungarian distance may be ideal for the Correspondence mechanism, but it is not defined for sets with different number of points, as it the case for the Re-sampling mechanism.

###### List of Tables

  1. 5.1 Summarization of the differences between datasets.
  2. 6.1 Mean norm of the flow vector of each dataset in the benchmark. A hypothetical model with zEPE = 1 on all datasets performs as well on Single ShapeNet as on Kitti, even though the EPE on those datasets is rather different. In a scene in KITTI most points belong to the ground, the mean norm is skwed to the average ego motion of the car. When the ground is removed, the mean norm is closer to the average motion of the dynamic objects in the scene.
  3. 6.2 Comparison of different methods performing flow estimation on KITTI. The best metrics of self-supervised methods are reported in bold. The underlined metrics indicate the best overall performance, regardless of the training method.
  4. 6.3 Performance of the baselines on the synthetic datasets.
  5. 6.4 FlowNet3D pretrained on FlyingThings3D, finetuned and evaluated on KITTI.
  6. 6.5 Flow Extractor trained with Adversarial Metric Learning and Cycle Consistency. Bold number show the best performing model per dataset.
  7. 6.6 Flow Extractor trained with Adversarial Metric Learning and Cycle Consistency. Models trained on Lyft and tested on both versions of KITTI.
  8. 6.7 Ablation studies on FlowNet3D trained with different auxiliary cycle consistency losses. The impact of different losses is quite expressive.
  9. 6.8 FlowNet3D trained on Adversarial Metric using different scaling factors for the Multi-Scale Triplet loss. The scaling factors are functions that use the level $l$ of the activations. The activations are organized from last to first. Zero only used the last activation.
  10. 6.9 Flow Extractors trained with nearest neighbor loss and cycle consistency loss.
  11. 6.10 EPE comparison on the impact of the mechanism used for training in the evaluation performance, FlyingThings3D used as dataset. All results are on supervised training, except for Point-PWC (self-sup) where we use the same self-supervised losses as the authors. * taken from [11], ** taken from [5].
  12. B.1 The architecture used for the Cloud Embedder. The output of every module was used on the Multi-Scale Triplet loss. Pointwise features are converted to global feature vectors by max pooling.
  13. B.2 The architecture used for the Cloud Embedder. The output of every module was used on the Multi-Scale Triplet loss. Point wise features are converted to global feature vectors by mean pooling.
  14. C.1 Difference in performance between Multi-Scale and shallow losses. The metrics are all better for the Multi-Scale loss. Even though the Multi-Scale loss adds some overhead to the calculations, the wall time is still lower.
  15. C.2 Difference in performance between Multi-Scale and shallow losses. Multi-Scale loss makes the training much more robust to reduced capacity of the Cluster Net.
  16. C.3 Difference in performance between Multi-Scale and shallow losses. Multi-Scale loss makes the training much more robust to reduced capacity of the Cluster Net.

### List of Algorithms

loa

##  Appendix A Algorithms

Input: index

// Load point cloud of an random object of ShapeNet $nframe=2$ ;

// number of frames

$npoints=1024$ ;

// number of points on output

$obj$ = ShapeNet.load(index) ;

// $obj.shape=[3,npoints]$

// Initialize parameters for initial transform ;

$\theta_{0}=[U(-\pi,\pi),U(-\pi,\pi),U(-\pi,\pi)]^{T}$ ;

// initial angle

$\mu_{0}=[U(5,6),U(5,6),U(5,6)]^{T}$ ;

// initial stretch factor

$T_{0}=[U(-1,1),U(-1,1),U(-1,1)]^{T}$ ;

// initial position

$M_{0}=\text{MakeRotationMatrix}(\theta_{0},\mu_{0})$

$obj_{0}=M_{0}\cdot obj+T_{0}$

// Initialize parameters for frame transform ;

$\theta=[U(-\frac{\pi}{15},\frac{\pi}{15}),U(-\frac{\pi}{15},\frac{\pi}{15}),U(-\frac{\pi}{15},\frac{\pi}{15})]^{T}$

$\mu=[U(0.9,1.1),U(0.9,1.1),U(0.9,1.1)]^{T}$

$T=[U(-0.25,0.25),U(-0.25,0.25),U(-0.25,0.25)]^{T}$

$M=\text{MakeRotationMatrix}(\theta,\mu)$

$t=1$

while  _$t\leq nframes$_ do

$obj_{t}=M\cdot obj_{t-1}+T$

$flow_{t-1}=obj_{t-1}-obj_{t}$

Draw $npoints$ random indices $idx$ without replacement from $obj_{t-1}$

$pc_{t-1}=obj_{t-1}[idx]$ ;

// select a random sample

$f_{t-1}=flow_{t-1}[idx]$ ;

// flow vectors of $pc_{t-1}$

end while

Output: $pc,f$

Algorithm 1 Single ShapeNet

There are minor implementation and computational differences between the
algorithms 1 and 2. Those modifications are in place to allow for the batch
computation of matrix transformations.

Input: index

// Load point cloud of an random object of ShapeNet $nframe=2$ ;

// number of frames

$npoints=1024$ ;

// number of points on output

$nobjs=U(2,20)$ ;

// number of objects

$n=1$

while  _$n\leq nobjs$_ do

$idx$ is randomly drawn; $obj_{n}$ = ShapeNet.load(idx) ;

$n+=1$

end while

// $obj.shape=[nbojs,npoints,3]$

// Initialize parameters for initial transform ;

$\theta_{0}=[U(-\pi,\pi),U(-\pi,\pi),U(-\pi,\pi)]^{T}$ ;

// initial angle

$\mu_{0}=[U(3,8),U(3,8),U(3,8)]^{T}$ ;

// initial stretch factor

$T_{0}=[U(-10,10),U(-10,10),U(-10,10)]^{T}$ ;

// initial position

// $theta_{0}.shape,\mu_{0}.shape,T_{0}.shape=[nobjs,3]$
$M_{0}=\text{MakeRotationMatrix}(\theta_{0},\mu_{0})^{T}$

// $M_{0}.shape=[nobjs,3,3]$ $obj_{0}=obj\cdot M_{0}+T_{0}$

// Initialize parameters for frame transform ;

$\theta=[U(-\frac{\pi}{15},\frac{\pi}{15}),U(-\frac{\pi}{15},\frac{\pi}{15}),U(-\frac{\pi}{15},\frac{\pi}{15})]^{T}$

$\mu=[U(0.9,1.1),U(0.9,1.1),U(0.9,1.1)]^{T}$

$T=[U(-0.25,0.25),U(-0.25,0.25),U(-0.25,0.25)]^{T}$

$M=\text{MakeRotationMatrix}(\theta,\mu)^{T}$

$t=1$

while  _$t\leq nframes$_ do

$obj_{t}=obj_{t-1}\cdot M+T$

$flow_{t-1}=obj_{t-1}-obj_{t}$

Draw $npoints$ random indices $idx$ without replacement from $obj_{t-1}$

$pc_{t-1}=obj_{t-1}[idx]$ ;

// select a random sample

$f_{t-1}=flow_{t-1}[idx]$ ;

// flow vectors of $pc_{t-1}$

$t+=1$

end while

Output: $pc,f$

Algorithm 2 Multi ShapeNet

The original FlyingThings3D [14] was designed to emulate stereo cameras. The
authors use a rendering engine to render diverse objects with independent
movement in a scene. Each scene was made of ten frames. For each frame the
following data is available:

  * •

RGB image of left and right camera;

  * •

Disparity map, which can be converted to a depth map via the focal length that
is provided;

  * •

Forward and backward optical flow;

  * •

Forward and backward disparity change;

The authors of [5] post-processed the original dataset into point clouds. We
use their post-processed dataset. The authors of [5] did not publish the
algorithm used for the post-processing, however, code was made available. We
summarize it in the algorithm 3.

Input: index

$npoints=1024$ ;

// number of points on output

$focalLength=1050$

$imgWidth,imgHeight=540,960$

Load $img_{0}$, $img_{1}$, $disparityMap_{0}$, $disparityMap_{1}$,
$disparityChange_{0}$, $opticalFlow_{0}$ from files

// The data is in grid format $[imgWidth,imgHeight,C]$

// $C$ is the number of channels (3 for RGB, 1 for disparity, 2 for optical
flow)

$depthMap_{0}=\frac{focalLength}{disparityMap_{0}}$

$xyGrid=\text{meshgrid}(\text{range}(imgWidth),\text{range}(imgHeight)$

$xyzMap_{0}[x]=(xyGrid[x]-(\frac{imgWidth}{2}))*\frac{focalLength}{depthMap_{0}}$

$xyzMap_{0}[y]=(xyGrid[y]-(\frac{imgHeight}{2}))*\frac{focalLength}{depthMap_{0}}$

$xyzMap_{0}[z]=depthMap_{0}$

$depthChangeMap=\frac{focalLength}{disparityMap_{0}}$

$xyNextGrid=xyGrid+opticalFlow_{0}$

$xyzChangeMap[x]=(xyNextGrid[x]-(\frac{imgWidth}{2}))*\frac{focalLength}{depthChangeMap}$

$xyzChangeMap[y]=(xyGrid[y]-(\frac{imgHeight}{2}))*\frac{focalLength}{depthChangeMap}$

$xyzChangeMap[z]=depthChangeMap$

$flowMap_{0}=xyzChangeMap-xyzChangeMap$

// Draw points for the first frame Draw random indices $idx$ from $img_{0}$

$rgb_{0}=img_{0}[idx]$

$pos_{0}=xyzMap_{0}[idx]$

$flow_{0}=flowMap_{0}[idx]$

$depthMap_{1}=\frac{focalLength}{disparityMap_{1}}$

$xyzMap_{1}[x]=(xyGrid[x]-(\frac{imgWidth}{2}))*\frac{focalLength}{depthMap_{1}}$

$xyzMap_{1}[y]=(xyGrid[y]-(\frac{imgHeight}{2}))*\frac{focalLength}{depthMap_{1}}$

$xyzMap_{1}[z]=depthMap_{1}$

Draw random indices $idx$ from $img_{1}$

$rgb_{1}=img_{1}[idx]$

$pos_{1}=xyzMap_{1}[idx]$

Output: $pos_{0},pos_{1},rgb_{0},rgb_{1},flow_{0}$

Algorithm 3 Point Cloud FlyingThings3D

##  Appendix B Experimental Setup

In this section we describe the experimental setup of individual experiments.

Tables 6.2 and 6.10: we used the same setup as made publicly available by the
authors of [11, 12] on [30, 31].

Table 6.3: non-leaning baselines have no training stage. The KNN estimator
uses $k=1$. We set $\gamma_{\operatorname*{\mathsf{knn}}}=1.0$ and
$\gamma_{\operatorname*{\mathsf{cc}}}=1.0$.

Table 6.3, supervised baselines:

  * •

Training Type: Supervised;

  * •

Flow Extractor training details:

    * –

Optimizer: Adam;

    * –

Learning Rate: $lr=0.0001$;

    * –

Betas: $\beta_{1}=0.0,\beta_{2}=0.99$;

    * –

Schedule: $0.5*lr$ every 5 epochs the evalaluation EPE does not decrease;

    * –

Weight Decay: $4\cdot 10^{-4}$;

    * –

Loss: l2.

  * •

Batch Size: $64$ for ShapeNet and Multi ShapeNet and $32$ for FlyingThings3D;

  * •

Stopping criteria: Training for at least 50 epochs, train concludes if no
improvement on the EPE is seen for consecutive 20 epochs or 500 epochs is
reached.

Table 6.4: finetunning perfomed with Kitti, used the following setup:

  * •

Training Type: Supervised, warm start;

  * •

Flow Extractor pretrained on FlyingThings3D. Training details:

    * –

Optimizer: Adam;

    * –

Learning Rate: $lr=0.0001$;

    * –

Betas: $\beta_{1}=0.0,\beta_{2}=0.99$;

    * –

Schedule: $0.5*lr$ every 5 epochs the evalaluation EPE does not decrease;

    * –

Weight Decay: $4\cdot 10^{-4}$;

    * –

Loss: l2.

  * •

Batch Size: $16$;

  * •

Stopping criteria: Training for at least 50 epochs, train concludes if no
improvement on the EPE is seen for consecutive 20 epochs or 500 epochs is
reached.

Table 6.9:

  * •

Training Type: Self-Supervised via KNN;

  * •

Flow Extractor training details:

    * –

Optimizer: Adam;

    * –

Learning Rate: $lr=0.0001$;

    * –

Betas: $\beta_{1}=0.0,\beta_{2}=0.99$;

    * –

Schedule: $0.5*lr$ every 5 epochs the evalaluation EPE does not decrease;

    * –

Weight Decay: $4\cdot 10^{-4}$;

    * –

Loss: $\mathcal{L}=\mathcal{L}_{knn}+\mathcal{L}_{cc}$

  * •

Batch Size: Batch Size: $64$ for ShapeNet and Multi ShapeNet and $32$ for
FlyingThings3D;

  * •

Stopping criteria: Training for at least 50 epochs, train concludes if no
improvement on the EPE is seen for consecutive 20 epochs or 500 epochs is
reached.

Tables 6.5, 6.7, 6.8:

  * •

Training Type: Self-Supervised via Adversarial Metric Learning;

  * •

Flow Extractor training details:

    * –

Optimizer: Adam;

    * –

Learning Rate: $lr=0.0005$;

    * –

Betas: $\beta_{1}=0.0,\beta_{2}=0.99$;

    * –

Schedule: $0.75*lr$ every 10 epochs the evalaluation EPE does not decrease;

    * –

Weight Decay: $4\cdot 10^{-4}$;

    * –

Loss: $\mathcal{L}=\mathcal{L}_{FE}$, defined by equation 4.8;

    * –

Cycle Consistency uses cosine similarity and l2 distance unless explicitly
stated.

  * •

Cloud Embedder training details:

    * –

FlowNet3D and PointPWC-net were trained with PointNet++ Cloud Embedder
described on Table B.2. The Segmenter was trained with PointNet Cloud Embedder
described on Table B.1.

    * –

Optimizer: Adam;

    * –

Learning Rate: $lr=0.00005$;

    * –

Betas: $\beta_{1}=0.0,\beta_{2}=0.99$;

    * –

Schedule: $0.75*lr$ every 10 epochs the evalaluation EPE does not decrease;

    * –

Weight Decay: $4\cdot 10^{-4}$;

    * –

Loss: $\mathcal{L}=\mathcal{L}_{CE}$, defined by equation 4.7.

  * •

Batch Size: Batch Size: $64$ for ShapeNet and Multi ShapeNet and $32$ for
FlyingThings3D;

  * •

Stopping criteria: Training for at least 50 epochs, train concludes if no
improvement on the EPE is seen for consecutive 40 epochs or 500 epochs is
reached.

PointNet Cloud Embedder
---
Module | Details
PointNet Feature Extractor 1 |  | layers: [64, 128, 256]
---
activation: ReLU
normalization: BatchNorm
output: point wise features
PointNet Feature Extractor 2 |  | layers: [1024, 1024, 512]
---
activation: ReLU
normalization: BatchNorm
output: global feature
MLP 1 |  | layers: [1024, 2048, 2048]
---
activation: ReLU
normalization: BatchNorm
output: global feature
MLP 2 |  | layers: [2048, 4096]
---
activation: ReLU
normalization: BatchNorm
output: global feature
Table B.1: The architecture used for the Cloud Embedder. The output of every
module was used on the Multi-Scale Triplet loss. Pointwise features are
converted to global feature vectors by max pooling. PointNet++ Cloud Embedder
---
Module | Details
PointNet++ Set Abstraction 1 |  | number of samples: 16
---
layers: [32, 64, 64]
activation: ReLU
normalization: BatchNorm
output: point wise features
PointNet++ Set Abstraction 2 |  | number of samples: 16
---
layers: [64, 128, 256]
activation: ReLU
normalization: BatchNorm
output: point wise features
PointNet++ Set Abstraction 3 |  | number of samples: 8
---
layers: [256, 512, 1024]
activation: ReLU
normalization: BatchNorm
output: point wise features
MLP 1 |  | layers: [128, 256]
---
activation: ReLU
normalization: BatchNorm
output: global feature
MLP 2 |  | layers: [512, 1024]
---
activation: ReLU
normalization: BatchNorm
output: global feature
MLP3 |  | layers: [2048, 4096]
---
activation: ReLU
normalization: BatchNorm
output: global feature
Table B.2: The architecture used for the Cloud Embedder. The output of every
module was used on the Multi-Scale Triplet loss. Point wise features are
converted to global feature vectors by mean pooling.

##  Appendix C Multi-Scale Triplet Loss

The first promising results came from experiments for flow extraction on point
clouds. The Multi-Scale Triplet Loss helps with the stability of the training
and boosts performance. The table C.1 gives a taste of those experiments. We
also performed classification and clustering experiments.

Loss type | EPE | Acc 01 | Acc 005 | wall time | epochs
---|---|---|---|---|---
Shallow Triplet L2 | 0.1622 | 26.57% | 6.37% | 4h | 250
Multi-Scale Triplet L2 | 0.1574 | 32.09% | 7.17% | 2h50min | 179
Table C.1: Difference in performance between Multi-Scale and shallow losses.
The metrics are all better for the Multi-Scale loss. Even though the Multi-
Scale loss adds some overhead to the calculations, the wall time is still
lower.

Setup: One network learns to map MNIST digits to a 2D space. We call it the
Cluster Network, it is trained using Hard Triplet Mining. The classification
is done via a linear Classifier, the classifier maps the 2D vector to a one-
hot encoded vector with 10 classes. The two models are trained together. For
evaluation it is considered:

  * •

Triplet Accuracy: the positive samples are mapped closer to the anchor than to
randomly chosen negative samples.

  * •

Classifier Accuracy: The Classifier trains for one epoch with the Cluster Net
frozen. I measure the overall accuracy of the classifier.

  * •

Time to 50% Classifier Accuracy: The number of epochs needed to reach 50% of
classification accuracy, that is to proxy how the first stage of training
evolves.

The motivation behind using a linear classifier is as follows. The Cluster Net
will ideally map the digits to be linearly separable. If that is the case, the
Classifier should get near 100% accuracy. However, if the Cluster Net is not
distinctive enough, it will make overlapping clusters. The Classifier Accuracy
will then drop.

The first experiment compares the same Cluster Network trained using deep and
shallow loss at difference capacities. The Network has 7 layers, 4 conv and 3
fully connected. The capacity is defined as the layer with the most number of
hidden units, all others are proportionally adjusted, a minimum of 4 hidden
units per layer is ensured.

Table C.2 shows how using Multi-Scale triplet loss helps to make training
robust to networks with different capacity.

Loss type | Capacity | Time to 50% | Classifier acc | Triplet Acc
---|---|---|---|---
Shallow | 512 | 35 | 73% | 96%
Multi-Scale (6) | 512 | 20 | 75% | 96%
Shallow | 32 | 30 | 71% | 93%
Multi-Scale (6) | 32 | 25 | 70% | 95%
Shallow | 4 | – | 44% | 86%
Multi-Scale (6) | 4 | 45 | 64% | 90%
Table C.2: Difference in performance between Multi-Scale and shallow losses.
Multi-Scale loss makes the training much more robust to reduced capacity of
the Cluster Net.

Another experiment shows the impact of the number of layers used for Multi-
Scale triplet loss. The number in the brackets indicates the number of layers
receiving triplet feedback. 6 means every layer receives a feedback, 3 means
every other layer and so on. For networks with limited capacity, dense
feedback seems useful too keep training progressing. If the network has enough
capacity, then it seems to be beneficial to allow for a few not-linear
transformations in between feedback points.

Loss type | Capacity | Time to 50% | Classifier acc | Triplet Acc
---|---|---|---|---
Shallow | 512 | 35 | 73% | 96%
Shallow | 4 | – | 44% | 86%
Multi-Scale (6) | 512 | 20 | 75% | 96%
Multi-Scale (6) | 4 | 45 | 64% | 90%
Multi-Scale (3) | 512 | 30 | 75% | 96%
Multi-Scale (3) | 4 | 50 | 61% | 89%
Multi-Scale (2) | 512 | 40 | 82% | 96%
Multi-Scale (2) | 4 | – | 46% | 84%
Multi-Scale (1) | 512 | 30 | 79% | 96%
Multi-Scale (1) | 4 | – | 44% | 83%
Table C.3: Difference in performance between Multi-Scale and shallow losses.
Multi-Scale loss makes the training much more robust to reduced capacity of
the Cluster Net.


# Flow-induced flapping of an inverted flag with non-uniform stiffness
distribution

Andres Goza1<EMAIL_ADDRESS>1Department of Aerospace Engineering,
University of Illinois at Urbana-Champaign, IL, USA

###### Abstract

We perform high-fidelity, two-dimensional (2D), fluid-structure interaction
(FSI) simulations at a Reynolds number of $Re=200$ of uniform flow past an
inverted flag (i.e., clamped at its trailing edge). The inverted flag system
can exhibit large-amplitude flapping motions (on the order of the flag length)
that can be converted to electricity via, e.g., piezoelectric materials. We
investigate the effect of structural nonuniformity in altering the FSI
dynamics compared with the uniform-stiffness scenario that has been thoroughly
characterized. We consider linear, quadratic, and cubic stiffness
distributions, and demonstrate that the FSI dynamics mirror those of a
uniform-stiffness flag with an appropriately defined effective stiffness. We
show that this effective stiffness can be computed simply via analysis of an
_in-vacuo_ Euler-Bernoulli beam. When expressed in terms of the effective
stiffness, the FSI dynamics of the nonuniform-stiffness flag exhibit the same
regimes—with many similarities in the detailed dynamics—as a uniform-stiffness
flag. This study opens questions about (i) what the optimal stiffness
distribution is for, e.g., energy harvesting capacity, and (ii) how to use
nonuniform (and possibly time-varying) stiffness distributions to control the
flag dynamics towards a desired state.

##  1 Introduction

Uniform flow past a deformable flag is commonly studied as an archetypal
fluid-structure interaction (FSI) system for a range of phenomena that are
pervasive to natural and engineered systems. Studies of the inverted-flag
system, where the flag is clamped at its trailing edge with respect to the
oncoming flow, have revealed that the inverted-flag system exhibits rich
dynamics with multiple behavioral regimes. One of these regimes involves
large-amplitude flapping (on the order of the flag’s length), which is
associated with large bending strains that provide potential for renewable
energy harvesting by, _e.g._ , coupling the flag to a piezoelectric material
(Shoele & Mittal, 2016; Orrego et al., 2017).

The parameters that govern the dynamics of the flag system are the
dimensionless mass and stiffness as well as the Reynolds number, defined here
(respectively) as

$M=\frac{\rho_{s}h}{\rho_{f}L},\;K=\frac{D}{\rho_{f}U^{2}L^{3}},\;Re=\frac{\rho_{f}UL}{\mu},$
(1)

where $\rho_{f}$ ($\rho_{s}$) is the fluid (structure) density, $U$ is the
freestream velocity, $L$ is the flag length, $\mu$ is the shear viscosity of
the fluid, $h$ is the flag thickness, and $D$ is the flexural rigidity of the
flag. A variety of behavioral regimes can be observed by varying these
dimensionless parameters (with the dimensionless stiffness $K$ playing the
most significant role in determining the dynamics).

For large stiffness values, an undeformed equilibrium state with the flag
remaining undeflected has been found experimentally, computationally, and
theoretically. This equilibrium has been shown to become unstable with
decreasing stiffness due to a divergence instability (Kim et al., 2013;
Gurugubelli & Jaiman, 2015; Sader et al., 2016a). This instability gives way
to a small-deflection stable state for low Reynolds numbers, $O(100-1{,}000)$
(Ryu et al., 2015; Gurugubelli & Jaiman, 2015). This deformed stable state was
shown to be a formal equilibrium of the full FSI system, and was found to give
way with a further decrease in stiffness to small-deflection flapping caused
by a supercritical Hopf bifurcation (Goza et al., 2018). These small-
deflection regimes have not been reported in higher Reynolds numbers
experiments, though this may be due in part to the narrow range of stiffness
values over which this regime is observed.

With further decrease in stiffness, large-amplitude flapping occurs with
flapping amplitudes that are commensurate with the flag length. This regime
has been linked for low mass ratios $M<1$ to a vortex-induced vibration (VIV;
see, _e.g._ , Sarpkaya (2004); Williamson & Govardhan (2004)), with flapping
initiated by leading-edge separation (Gurugubelli & Jaiman, 2015) and
subsequent synchronization occurring between the integrated fluid forces and
flapping dynamics occurring near a Strouhal number of $0.2$ (Sader et al.,
2016a). Disruptions to this VIV behavior have been observed. Experimental and
numerical investigations have identified that even when serrations and
splitter plates were used to disrupt key flow structures, flapping persisted
at different frequencies and produced distinct vortex shedding processes
(Tavallaeinejad et al., 2020b; Gurugubelli & Jaiman, 2019). Moreover, for
large mass ratios ($M>1$) flapping is qualitatively distinct from the
classical VIV phenomena, with simulations highlighting different flapping-
vortex shedding interactions (Goza et al., 2018) and theoretical work
demonstrating that a quasi-steady fluid model—devoid of vortex
shedding—predicted many features of the heavy-flag flapping dynamics
(Tavallaeinejad et al., 2020a).

Further decrease in stiffness yields to a de-synchronization between flapping
and vortex shedding (Goza et al., 2018) and subsequent chaotic flapping, which
has been found experimentally at high Reynolds numbers of $O(10^{4})$ (Sader
et al., 2016a) and computationally at the lower Reynolds number of $200$ (Goza
et al., 2018). Goza et al. (2018) characterized this regime as one with an
attractor involving (seemingly) random switching between large-amplitude
flapping and the deformed mode. The authors demonstrated that chaotic flapping
does not occur at low Reynolds numbers $\lesssim 50$ (where classical bluff-
body vortex shedding is absent) or for massive flags ($M>1$), which suggests
that chaos arises from a de-synchronization/breakdown of VIV behavior in
large-amplitude flapping of light flags, but before the deformed mode becomes
a global attractor.

A further decrease in stiffness results in the deformed mode regime, in which
there is small-amplitude flapping occurs about a large mean-deflection state.
This regime has been observed experimentally (Kim et al., 2013) and
computationally (Ryu et al., 2015; Gurugubelli & Jaiman, 2015; Shoele &
Mittal, 2016). For Reynolds numbers below $\lesssim 50$ (where classical
bluff-body vortex shedding ceases), the regime is an equilibrium state of the
FSI system (Goza et al., 2018). At higher Reynolds numbers, flapping occurs
near the classical bluff-body frequency scaling of 0.2 (Shoele & Mittal,
2016).

The sensitivity of the large-amplitude flapping regime highlights a challenge
for energy harvesting: small changes in dimensionless stiffness ($K$,
achievable through, _e.g._ , a change in the oncoming flow velocity) can
elicit a transition from large-amplitude flapping to a different regime that
is less useful for energy capture. A possibility for addressing this
sensitivity is to manipulate the flag’s structural properties (made possible
by advances in modern materials) to alter the flag dynamics. To inform how to
tune these material properties, a key question is what role structural
nonuniformity plays in dictating the dynamics of the inverted-flag system. We
will answer this question here in the case of a flag with a linear, quadratic,
and cubic stiffness distribution. Our studies use high-fidelity two-
dimensional (2D) simulations at $Re=200$, $M=0.5$. We demonstrate that the
dynamics of the nonuniform-stiffness flag can be predicted by a uniform-
stiffness flag, provided that an effective stiffness is utilized. We further
show that this effective stiffness is different from the mean stiffness, and
can be predicted by a simple calculation involving an Euler-Bernoulli beam in
a vacuum. When expressed in terms of the effective stiffness, the FSI dynamics
of the nonuniform-stiffness flag are shown to exhibit the same regimes—with
many similarities in the detailed dynamics—as a uniform-stiffness flag.

We restrict our attention to the aforementioned mass ratio $M$ and Reynolds
number $Re$. The relative robustness of the dynamics to these two parameters
suggests that features of our conclusions will hold for other flag masses and
Reynolds numbers. We also do not investigate 3D effects here. The commonality
between 3D experiments and 2D simulations highlighted above as well as
systematic studies of aspect ratio suggest that the majority of 2D phenomena
persist in 3D, provided that the aspect ratio is not too small (Sader et al.,
2016b; Tavallaeinejad et al., 2021). Similarly, the inclination angle of the
flow with respect to the flag is not considered here, but can alter the flag
dynamics considerably (Shoele & Mittal, 2016; Huertas-Cerdeira et al., 2021).
An investigation into the role of 3D phenomena and the initial inclination
angle is left to future study.

##  2 Problem statement and numerical method

In this work, we consider an inverted-flag system in which the flag is clamped
at its leading edge with respect to the oncoming flow. The dimensionless mass
and Reynolds number are fixed in this study at $M=0.5$ and $Re=200$,
respectively. Drawing inspiration from Floryan & Rowley (2020), who studied
the role of structural nonuniformity in bio-inspired locomotion, we express
the dimensionless stiffness distribution in terms of Legendre polynomials as
$K(s)=2\overline{K}+\sum_{j=1}^{n}c_{j}L_{j}(s)$, where $s\in[-1,1]$ is a
parametrization of the flag (non-dimensionalized by the flag length $L$),
$\overline{K}$ is the mean stiffness, and $L_{j}$ is the $j^{th}$ Legendre
polynomial. This stiffness distribution uses the same non-dimensionalization
as in (1). In this work, we consider linear, quadratic, and cubic stiffness
distributions, corresponding to $n=1,$ $2$, and $3$, respectively. The benefit
of the representation via Legendre polynomials is that $\langle
L_{j},1\rangle=0$, where $\langle\cdot,\cdot\rangle$ is the standard
(unweighted) inner product on the interval $[-1,1]$. That is, varying the
coefficients $c_{1}$, $c_{2}$, and $c_{3}$ does not affect the mean stiffness
$\overline{K}$. This approach allows for the effect of stiffness distribution
to be investigated without inheriting effects from changing the mean
stiffness. The only restriction on $c_{1}$, $c_{2}$, and $c_{3}$ is thus that
$K(s)>0$ $\forall s\in[-1,1]$.

To meaningfully sample the parameter space associated with the various
coefficients, we focus on mean stiffness values of $\overline{K}=0.15$,
$0.25$, $0.35$, and $0.45$. These mean stiffnesses correspond in the uniform-
distribution studies to the deformed mode, chaotic, large-amplitude, and
small-deflection flapping regimes, respectively, for this choice of mass ratio
and Reynolds number. For each mean stiffness $\overline{K}$ and polynomial
degree $n$, the coefficients $c_{1}$, $c_{2}$, and $c_{3}$ are varied as
follows. We first construct the feasible set of coefficients for the given
mean stiffness $\overline{K}$ and polynomial degree $n$. We then perform
simulations using a uniformly distributed set of coefficients within the
feasible set. Figure 1 shows the feasible set and values of $c_{1}$, $c_{2}$
considered for $\overline{K}=0.35$ and $n=2$. The specific number of
simulations varied between fifteen and forty for each choice of $\overline{K}$
and $n$ based on the size of the feasible set. In total, seventy four, eighty
three, and eighty nine runs were performed for $n=1$, $2$, and $3$,
respectively.

Figure 1: Example for $\overline{K}=0.35$, $n=2$ of how the coefficients
$c_{1}$, $c_{2}$, and $c_{3}$ were selected to define the nonuniform stiffness
distribution $K(s)$. Left subplot: feasible set of coefficients (shaded
region) and coefficients that were used for the simulations (red markers).
Right subplot: stiffness distributions for a subset of the utilized
coefficients (green open circles in the left subplot). The coefficients used
for other values of mean stiffness and polynomial order were obtained through
the same process.

All simulations are run using the immersed boundary method of Goza & Colonius
(2017). The method has been validated on a number of flapping flag problems
for flags in both the conventional and inverted configuration. The method was
also used to identify physical mechanisms in the inverted-flag system for a
flag of uniform stiffness (Goza et al., 2018). The simulation parameters for
the present study are the same as in the cited reference; a convergence test
therein demonstrates the suitability of those simulation parameters (a
separate convergence test was run using the parameters for the present study
and found to give similar results—this outcome is expected, as the dynamics
for this non-uniform stiffness setting will be shown below to be commensurate
with those from the uniform-stiffness case).

For all simulations, the flag is initialized in the undeformed configuration,
and dynamics are triggered by introducing a small body force at early time. To
ensure that long-term behavior is captured, the simulations are run for a
minimum of 150 convective time units (for the cases exhibiting chaotic
dynamics, simulations are run for 300 convective time units). In the cases
exhibiting limit-cycle dynamics, this final time amounts to several dozen
periods being captured; for reference, in the left two columns of figure 2
below, the $x$-limits are cut off at 100 convective time units.

##  3 Summary for a uniform flag and intuition behind effective stiffness

To build context for the ensuing non-uniform stiffness results, we review here
the salient dynamics of a uniform-stiffness inverted flag. Figure 2
demonstrates the various regimes observed by decreasing the stiffness $K$. The
frequency content (bottom right subplot) highlights that the transition to
large-amplitude flapping coincides (at this mass ratio of $M=0.5$) with a jump
from low flapping frequency to one that is synchronized to vortex shedding at
the classical bluff body scaling of 0.2. This is one of the markers of the VIV
behavior that the system exhibits (Sader et al., 2016a). With decreasing
stiffness, the flag frequency decreases steadily, signaling a de-
synchronization between flapping and vortex shedding. Eventually, the de-
synchronization disrupts large-amplitude flapping and triggers chaotic
flapping with broadband frequency content. Finally, with further decreases in
stiffness, the flag enters into the deformed mode regime, and the dominant
frequency is near 0.2 and its harmonics (though there are relatively small
peaks for some stiffnesses at lower frequencies, $\lesssim 0.04$). The
prominent frequency signatures of 0.2 reflect that the small-amplitude
flapping is driven by forcing from canonical bluff-body vortex shedding
(Shoele & Mittal, 2016). See Goza et al. (2018) for a more detailed
characterization of the uniform-stiffness dynamics for this choice of Reynolds
number, mass ratio, and collection of stiffness values.

Figure 2: Left two columns: transverse displacement of the flag’s leading
edge, $d_{y}$, as a function of time for the undeformed equilibrium ( ),
small-deflection equilibrium ( ), small-deflection flapping ( ), large-
amplitude flapping ( ), chaotic flapping ( ), and the deformed mode ( )
regimes. Regime transitions were triggered by decreasing the dimensionless
stiffness $K$. Third column: Alternative view of the regimes shown as a
bifurcation diagram (top row) and frequency content of the tip displacement
versus stiffness $K$ (bottom row). For the bifurcation diagrams, at a given
stiffness the markers show the leading-edge displacement at zero velocity.
(For example, for a stiffness in the chaotic regime, the markers at various
$d_{y}$ values reflect irregular direction changes; _c.f._ , the corresponding
timetrace plot in the second column). The frequency information in the bottom
right subplot provides, for each $K$, the frequency peaks from a power
spectral density analysis with at least 20% of the frequency peak with maximum
energy. To indicate the spread of energy about each frequency peak, the color
extends from the peak across neighboring frequencies with at least 5% of the
peak energy. (For example, for $K^{-1}\approx 5$ there are two dominant
frequencies, near $0.03$ and $0.1$, and there is considerable spread in energy
across neighboring frequencies for the peak near $0.03$).

We now provide archetypal results in figure 3 for a linear stiffness
distribution ($c_{2}=c_{3}=0$) with mean stiffness $\overline{K}=0.35$ and
various $c_{1}$ values. Note that since $L_{1}(x)=x$, $c_{1}$ is the slope,
with negative (positive) values corresponding to stiffness being more isolated
towards the leading (trailing) edge. The figure demonstrates that $c_{1}$ acts
as a bifurcation parameter, analogous to the role that $K$ plays for the
uniform stiffness case; _c.f._ , figure 2. Indeed, the linearly distributed
case of figure 3 shows that the same regimes (undeformed, small-deflection
equilibrium, small-deflection flapping, large-amplitude flapping, and deformed
mode) exist, and that the bifurcation sequence is the same with decreasing
$c_{1}$ as for decreasing $K$. Moreover, the frequency content of the linearly
distributed case (bottom right subplot of figure 3) demonstrates that the same
dynamics are present as in the uniform-stiffness case: small-amplitude
flapping of low frequency gives way to VIV flapping at large amplitudes, and a
subsequent de-synchronization of flapping and vortex shedding yields chaotic
flapping and eventually the deformed-mode regime characterized by some
broadband low frequency dynamics, but even stronger frequency peaks near $0.2$
and its integer harmonics.

Figure 3: Analog of figure 2 except $c_{1}$ is varied, not $K$.

These facts suggest that the flag with a linear stiffness distribution acts as
a uniform stiffness flag with an effective stiffness different from the mean
stiffness: for a given mean stiffness, flags with a stiffer trailing edge
($c_{1}>0$) act effectively stiffer and flags with a more flexible trailing
edge ($c_{1}<0$) act more flexible than a uniform stiffness flag with
stiffness $K\equiv\overline{K}$.

We will define this effective stiffness from an analogous uniform flag. To
build intuition, we first consider an Euler-Bernoulli beam in a vacuum with
mean stiffness $\overline{K}=0.35$ and linear stiffness distribution
($c_{2}=c_{3}=0$). This beam model is a linearized variant of the
geometrically nonlinear structural model used in our numerical method,
described in section 2. We show in figure 4 the leading eigenvectors and
eigenvalues (computed by an eigendecomposition of a prototypical finite
element discretization) for various values of $c_{1}$. The figure demonstrates
that while the second through sixth eigenvectors vary in shape with varying
$c_{1}$, to reasonable approximation the leading eigenvector does not. The
similarity of the leading eigenvector with varying $c_{1}$ suggests that the
dynamics can be represented by a uniform-stiffness structure, as the FSI
inverted-flag system is dominated by first mode dynamics. The right subplot
shows that the natural frequency increases (i.e., the effective stiffness
increases) as stiffness is focused more towards the clamped (trailing) edge.
This trend suggests a procedure for computing the effective stiffness, which
we employ here. Given a stiffness distribution, we define the effective
stiffness as the stiffness for a uniform-property flag that has the same
leading natural frequency as the non-uniform stiffness flag of interest. To
indicate why focusing stiffness towards the trailing edge (leading edge) leads
to an effectively stiffer (more flexible) flag, we note that the strain energy
within the flag is given by $\int_{-1}^{1}K(s)\kappa(s)ds$, where $\kappa(s)$
represents the curvature of the flag. The material stiffness $K(s)$ thus
contributes more to the strain energy—and thereby overall stiffness—in areas
where curvature is more pronounced. The clamped end is a natural location to
expect these regions of high curvature.

Figure 4: Leading eigenvectors (left two columns) and natural frequencies,
$\omega^{n}_{j}$, versus index $j$ (right column) for an _in-vacuo_ Euler-
Bernoulli beam with mean stiffness $\overline{K}=0.35$ and linear stiffness
distribution ($c_{2}=c_{3}=0$). The insert in the natural frequency plot shows
the natural frequency of the leading mode, index $j=1$.

We note that the interaction of the flag with the surrounding flow will alter
the natural frequency of the FSI system. We do not argue here that the vacuum-
scaled natural frequency is the correct natural frequency for the FSI system,
but instead that the effect the flow has on the non-uniform flag is similar to
its effect on a uniform flag when an appropriately defined effective stiffness
(distinct from the mean stiffness) is used.

##  4 Representing dynamics with the effective stiffness

We now utilize the definition of the effective stiffness, $K_{e}$, to
represent the dynamics for linear, quadratic, and cubic stiffness
distributions in figure 5. The left column of figure 5 is created from two
hundred forty six runs, using all three distributions and mean stiffness
values of $\overline{K}=0.15$, $0.25$, $0.35$, and $0.45$ (corresponding in
the uniform case to the deformed mode, chaotic, large-amplitude, and small-
deflection flapping regimes, respectively). The left column of the figure
demonstrates that the regime transitions occur at nearly the same $K_{e}$ as
for the uniform-stiffness case (c.f., figure 2). In addition, the frequency
content suggests that the transitions transpire by the same mechanisms with
decreasing effective stiffness—a jump in frequency to large-amplitude VIV
flapping, de-synchronization that triggers chaos and broadband frequency
dynamics, and the eventual settling of the flag to one side in a deformed-mode
state characterized by small-amplitude flapping due to flow forcing from
bluff-body vortex shedding.

The most notable exceptions to the predictive capabilities of $K_{e}$ are the
occurrence of chaotic dynamics where the deformed mode is predicted, and the
difference of flapping amplitude across the large-amplitude flapping regime.
For the former exception, the corresponding frequency content demonstrates
that there are prominent signatures at 0.2 and its harmonics, which reflects
the fact that the flag spends significant time in a quasi-deformed mode state
before flapping to the other side of the centerline.

Figure 5: Left column: analog of the right column of figure 2 for a variety of
linear, quadratic, and cubic non-uniform stiffness distributions. For ease of
comparison with the uniform-stiffness case, the color scheme is the same as in
figure 2. Right column: time trace of $d_{y}$ for five cases with an effective
stiffness of $K_{e}\approx 0.35$ (the legend provides the corresponding row in
figure 6, where parameter details are given).

Within the large-amplitude flapping regime, for all but one of the cases the
dynamics evolve onto a limit cycle with a fixed amplitude (to within $5\%$).
Thus, the variation in tip displacement $d_{y}$ is because different stiffness
distributions can lead to different flapping amplitudes—a notable contrast
from the nearly constant amplitude across this regime for the uniform case.
The periodic dynamics and variation in flapping amplitude are demonstrated
further in the right column of figure 5. Yet, despite this apparent difference
from the uniform case, there remain marked similarities in the detailed
dynamics. Figure 6 provides vorticity snapshots at equally spaced instances
during a flapping cycle for the five stiffness distributions from the
timetrace plot in the right column of figure 5. All five cases produce similar
vortex dynamics, with a leading-edge vortex (LEV) and trailing-edge vortex
(TEV) formed when the flag is at its peak and minimal displacement values,
respectively. Moreover, a P+P wake structure forms and advects in nearly
identical fashion. At the same time, for the instances of peak displacement
the case with smaller flapping amplitude leads to significant additional
bending near the chordwise center. This diversity of bending within common FSI
dynamics raises questions about what bending distributions are most useful for
energy capture within large-amplitude flapping.

Figure 6: Each row provides snapshots of vorticity, $\omega$, at four equally
spaced time instances across a flapping period for a flag with effective
stiffness near $K_{e}=0.35$. First row: $\overline{K}=0.35$,
$c_{1}=c_{2}=c_{3}=0$; second row: $\overline{K}=0.25$, $c_{1}=0.1765,$
$c_{2}=0.0588$, $c_{3}=0$, $K_{e}=0.350$; third row: $\overline{K}=0.45$,
$c_{1}=0.5172,$ $c_{2}=0.7241$, $c_{3}=0$, $K_{e}=0.351$; fourth row:
$\overline{K}=0.45$, $c_{1}=0.1765,$ $c_{2}=0.0588$, $c_{3}=-0.05882$,
$K_{e}=0.351$; fifth row: $\overline{K}=0.45$, $c_{1}=-0.1579,$
$c_{2}=0.1579$, $c_{3}=0.2632$, $K_{e}=0.351$. There are thirty four vorticity
contour levels over $\omega\in[-10,10].$

##  5 Conclusions

We performed high-fidelity 2D FSI simulations of flow past an inverted flag
with linear, quadratic, and cubic distributions. The stiffness distribution
was defined using Legendre polynomials so that stiffness distribution could be
systematically varied without affecting the mean stiffness. We demonstrated
that the dynamics of the non-uniform flag could be represented as those of a
uniform flag with an appropriately defined effective stiffness. This effective
stiffness was computed in a straightforward manner using an _in-vacuo_ Euler-
Bernoulli beam: the effective stiffness was defined as the stiffness of a
uniform flag with the same leading natural frequency as that of the non-
uniform stiffness flag. When cast in terms of this effective stiffness, the
same regimes (with minor exceptions) occurred for linear, quadratic, and cubic
distributions as in the uniform case. Moreover, many of the same dynamical
details persisted. Large-amplitude flapping continued to be characterized by
synchronization between vortex shedding and flapping, chaotic flapping was
triggered by a de-synchronization of these phenomena, and the deformed mode
persisted as a state in which the flag passively responded to bluff-body
driven vortex shedding. Moreover, within large-amplitude flapping the same
vortex-shedding and advection processes were found to persist even when the
flapping amplitude varied by up to roughly $25\%$ from the uniform-stiffness
flag. This outcome demonstrates that the bending profile under these different
distributions is modulated to accommodate consistent FSI dynamics across a
relatively wide range of stiffness distribution. This investigation opens
questions about which bending profiles within large-amplitude flapping are
optimal for energy capture, and how to incorporate a (potentially time
varying) stiffness distribution to obtain desired FSI dynamics. This study
considered a specific mass ratio and did not include an investigation into the
effect of initial inclination angle of the flag. A characterization of these
effects as well as the role of three dimensionality are avenues for future
study.

##  6 Declaration of interest

The author reports no conflict of interest.

## References

  * Floryan & Rowley (2020) Floryan, Daniel & Rowley, Clarence W 2020 Distributed flexibility in inertial swimmers. Journal of Fluid Mechanics 888.
  * Goza & Colonius (2017) Goza, Andres & Colonius, Tim 2017 A strongly-coupled immersed-boundary formulation for thin elastic structures. Journal of Computational Physics 336, 401–411.
  * Goza et al. (2018) Goza, Andres, Colonius, Tim & Sader, John E 2018 Global modes and nonlinear analysis of inverted-flag flapping. Journal of Fluid Mechanics 857, 312–344.
  * Gurugubelli & Jaiman (2015) Gurugubelli, PS & Jaiman, RK 2015 Self-induced flapping dynamics of a flexible inverted foil in a uniform flow. Journal of Fluid Mechanics 781, 657–694.
  * Gurugubelli & Jaiman (2019) Gurugubelli, Pardha S & Jaiman, Rajeev K 2019 Large amplitude flapping of an inverted elastic foil in uniform flow with spanwise periodicity. Journal of Fluids and Structures 90, 139–163.
  * Huertas-Cerdeira et al. (2021) Huertas-Cerdeira, Cecilia, Goza, Andres, Sader, John E, Colonius, Tim & Gharib, Morteza 2021 Dynamics of an inverted cantilever plate at moderate angle of attack. Journal of Fluid Mechanics 909.
  * Kim et al. (2013) Kim, Daegyoum, Cossé, Julia, Cerdeira, Cecilia Huertas & Gharib, Morteza 2013 Flapping dynamics of an inverted flag. Journal of Fluid Mechanics 736.
  * Orrego et al. (2017) Orrego, Santiago, Shoele, Kourosh, Ruas, Andre, Doran, Kyle, Caggiano, Brett, Mittal, Rajat & Kang, Sung Hoon 2017 Harvesting ambient wind energy with an inverted piezoelectric flag. Applied energy 194, 212–222.
  * Ryu et al. (2015) Ryu, Jaeha, Park, Sung Goon, Kim, Boyoung & Sung, Hyung Jin 2015 Flapping dynamics of an inverted flag in a uniform flow. Journal of Fluids and Structures 57.
  * Sader et al. (2016a) Sader, John E, Cossé, Julia, Kim, Daegyoum, Fan, Boyu & Gharib, Morteza 2016a Large-amplitude flapping of an inverted flag in a uniform steady flow–a vortex-induced vibration. Journal of Fluid Mechanics 793.
  * Sader et al. (2016b) Sader, John E, Huertas-Cerdeira, Cecilia & Gharib, Morteza 2016b Stability of slender inverted flags and rods in uniform steady flow. Journal of Fluid Mechanics 809, 873–894.
  * Sarpkaya (2004) Sarpkaya, Turgut 2004 A critical review of the intrinsic nature of vortex-induced vibrations. Journal of Fluids and Structures 19, 389–447.
  * Shoele & Mittal (2016) Shoele, Kourosh & Mittal, Rajat 2016 Energy harvesting by flow-induced flutter in a simple model of an inverted piezoelectric flag. Journal of Fluid Mechanics 790, 582–606.
  * Tavallaeinejad et al. (2020a) Tavallaeinejad, Mohammad, Païdoussis, Michael, Legrand, Mathias & Kheiri, Mojtaba 2020a Instability and the post-critical behaviour of two-dimensional inverted flags in axial flow. Journal of Fluid Mechanics 890.
  * Tavallaeinejad et al. (2020b) Tavallaeinejad, Mohammad, Païdoussis, Michael, Salinas, Manuel Flores, Legrand, Mathias, Kheiri, Mojtaba & Botez, Ruxandra 2020b Flapping of heavy inverted flags: a fluid-elastic instability. Journal of Fluid Mechanics 904.
  * Tavallaeinejad et al. (2021) Tavallaeinejad, Mohammad, Salinas, Manuel Flores, Païdoussis, Michael P, Legrand, Mathias, Kheiri, Mojtaba & Botez, Ruxandra M 2021 Dynamics of inverted flags: Experiments and comparison with theory. Journal of Fluids and Structures 101, 103199.
  * Williamson & Govardhan (2004) Williamson, Charles HK & Govardhan, R. 2004 Vortex-induced vibrations. Annual Review of Fluid Mechanics 36, 413–455.


11institutetext: UC San Diego, La Jolla, CA, USA 92093 22institutetext:
Princeton University, Princeton, NJ, USA 08544 33institutetext: Cornell
University, Ithaca, NY, USA 14853 44institutetext: Fermilab, Batavia, IL, USA
60510-5011 55institutetext: US Naval Academy, Annapolis, MD, USA 21402

# Generalizing mkFit and its Application to HL-LHC

Giuseppe Cerati 44 Peter Elmer 22 Patrick Gartung 44 Leonardo Giannini 11
Matti Kortelainen 44 Vyacheslav Krutelyov 11 Steven Lantz 33 Mario
Masciovecchio 11 Tres Reid 33 Allison Reinsvold Hall 55 Daniel Riley 33 MatevÅ
Tadel 11<EMAIL_ADDRESS>Emmanouil Vourliotis 11 Peter Wittich 33 Avi Yagil on
behalf of the CMS collaboration 11

###### Abstract

mkFit is an implementation of the Kalman filter-based track reconstruction
algorithm that exploits both thread- and data-level parallelism. In the past
few years the project transitioned from the R&D phase to deployment in the
Run-3 offline workflow of the CMS experiment. The CMS tracking performs a
series of iterations, targeting reconstruction of tracks of increasing
difficulty after removing hits associated to tracks found in previous
iterations. mkFit has been adopted for several of the tracking iterations,
which contribute to the majority of reconstructed tracks. When tested in the
standard conditions for production jobs, speedups in track pattern recognition
are on average of the order of 3.5x for the iterations where it is used (3-7x
depending on the iteration). Multiple factors contribute to the observed
speedups, including vectorization and a lightweight geometry description, as
well as improved memory management and single precision. Efficient
vectorization is achieved with both the icc and the gcc (default in CMSSW)
compilers and relies on a dedicated library for small matrix operations,
Matriplex, which has recently been released in a public repository. While the
mkFit geometry description already featured levels of abstraction from the
actual Phase-1 CMS tracker, several components of the implementations were
still tied to that specific geometry. We have further generalized the geometry
description and the configuration of the run-time parameters, in order to
enable support for the Phase-2 upgraded tracker geometry for the HL-LHC and
potentially other detector configurations. The implementation strategy and
high-level code changes required for the HL-LHC geometry are presented.
Speedups in track building from mkFit imply that track fitting becomes a
comparably time consuming step of the tracking chain. Prospects for an mkFit
implementation of the track fit are also discussed.

##  1 Introduction

The mkFit project was started in 2014 with the goal of exploring how the
traditional Kalman filter based track fitting and track finding kalmanfit  can
be rethought and optimized in the age of – at that time, novel – many-core,
vectorized computing architectures. After initial positive results on
simplified detector geometries the focus was shifted to applying the mkFit
algorithm to silicon detector track finding for the CMS experiment cms-det ,
resulting in a viable prototype in 2018 and culminating in a final one-year
integration and validation campaign in 2021. Since the beginning of LHC Run3
in 2022 cms-run3 , mkFit is used by CMS to reconstruct five out of twelve
tracking iterations covering 90% of found tracks with
$p_{T}>0.5\,\rm{G}e\rm{V}/c$. An in-depth review of mkFit, including detailed
motivation and algorithm description, has been published mkfit-ref . An
overview of early work with further references can be found in mkfit-chep18
and physics performance of mkFit in CMS Run3 is available as a CMS Detector
Performance note mkfit-dpnote-2022 .

This paper focuses on improvements and extensions of mkFit that were required
to support running of multiple tracking iterations in _CMS software_ (CMSSW)
as well as to prepare it for implementation of tracking after the Phase-2
detector upgrades expected around 2030, in the _High Luminosity LHC_ (HL-LHC)
era. Section 2 introduces how mkFit is structured and run within CMSSW. A
detailed presentation of required generalizations of geometry description,
configuration and steering systems is given in section 3. Currently ongoing
and planned or possible future work is discussed in section 4.

##  2 mkFit in CMSSW

mkFit was initially developed as a standalone tracking library and at first
included into CMSSW as an external package. This mode of operation was used
for development, physics performance tuning and benchmarking. However, one of
the conditions for using mkFit in production was for the code to be
incorporated into the core CMSSW distribution, to make the software building,
configuration, and patching for online and offline use compatible with CMS’s
requirements for computing operations. This section discusses the high-level
code structure of mkFit in CMSSW; outlines steps performed by mkFit in a
typical CMSSW reconstruction job; and, finally, presents some highlights of
physics and computing performance.

###  2.1 Code structure

mkFit code is structured into three CMSSW packages:

  * •

MkFitCore holds the central components of mkFit, including all computational
algorithms, internal data formats and geometry description, as well as
configuration structures and related processing code. This core package is
independent of any experiment or geometry details. It does not depend on or
interact directly with any CMSSW modules or data formats.

  * •

MkFitCMS contains helper algorithms, called _standard functions_ , that
perform specific tasks during track finding: seed pre-processing, candidate
scoring, candidate filtering, and duplicate removal. These codes use mkFit
internal data formats and still do not depend on CMSSW.

  * •

MkFit is the actual bridge between mkFit and CMSSW. It defines the CMSSW
producer modules for both configuration and data-processing. It uses CMSSW
specific mechanisms to pull in configuration and event-data, transforms them
into mkFit internal structures and calls appropriate steering functions. It
depends both on MkFitCore (data-formats and geometry) and MkFitCMS (steering
and standard functions) packages.

Standalone operation of mkFit is still possible and is frequently used for
validation, tuning, development, and debugging. To support this mode, packages
MkFitCore and MkFitCMS contain additional code and makefiles in sub-directory
standalone/ that is not used by CMSSW build or touched by CMSSW code
managements tools. This allows for keeping all mkFit related files stored in a
single repository. A minimal additional repository with external packages that
would otherwise be used from CMSSW or CMSSW’s external software still needs to
be maintained separately for standalone
builds.111https://github.com/trackreco/mkFit-external

###  2.2 Track finding algorithm

Details of CMS track reconstruction and iterative tracking can be found in
cms-tracking . Here we are concerned with processing as it occurs for every
iteration after the seed tracks have been found. Input to track finding is a
vector of seed tracks, each consisting of a list of associated hits and the
initial estimate of the track parameters at the final, outermost point. After
that, mkFit processing steps are as follows:

  1. 1.

_Seed cleaning_ is performed. As mkFit processes seeds in parallel it can not
rely on hit masking in order to exclude seeds whose hits have already been
consumed by previously found tracks.

  2. 2.

_Seed partitioning_ reshuffles the seeds into _tracking regions_ (barrel,
transition, and end-cap). Those define the sequence in which detector layers
will be visited. Additionally, the seeds are sorted in ${\eta,\varphi}$-space
to improve hit access coherency during later steps.

  3. 3.

_Forward search_ proceeds through the detector layers for the given tracking
region going outwards from the seeding layers. For each seed, combinatorial
search with a limit on the maximum concurrent number of candidates is
performed, adding new hits on each layer while allowing for a limited number
of missed layers and a single additional “detector overlap” hit. At the end,
either because the edge of the tracker is reached or because no more new hits
are found, the best-scoring candidate is chosen as the representative.
Optionally, a quality filter can be applied before the next step.

  4. 4.

_Backward fit_ re-traverses each found track backwards, refitting the track
parameters. If the seeding region for the current iteration does not extend
all the way to the vertex region, a combinatorial _backward search_ can also
be performed, going inwards from the first known hit.222Optionally, some or
all of the seed-region hits can be dropped and new hits searched for in the
seeding layers as well. Again, the search is stopped when the innermost layers
of the detector have been reached or if no new hits are found for a given
seed. If the search has been performed, the best candidate is chosen as the
final representative.

  5. 5.

_Quality filtering & duplicate removal_ are performed on the resulting tracks.

After the track finding for each iteration is complete, two more steps are
performed by other CMSSW modules: final fit including outlier rejection; and
final quality selection, based on a multivariate algorithm.

In standalone operation, the same mkFit processing steps are performed. Pre-
processed input seeds and hit data are read from a custom binary file. Final
fit and quality selection are not performed, but there is a standalone version
of validation comparing the found tracks to simulated ones.

###  2.3 Physics and computational performance

Here we present two highlights from the CMS Detector Performance note mkfit-
dpnote-2022 .333In this note mkFit has also been used for PixelLess iteration
that was later removed due to poor performance for low-momentum highly
displaced tracks. Both cases compare relevant quantities before and after the
inclusion of mkFit in the standard Run3 track reconstruction of simulated
$t\bar{t}$ events with an average pile-up of 65.

Figure 1: Comparison of physics performance of legacy track reconstruction
(red markers) and the new Run3 tracking configuration where mkFit is used for
6 of twelve tracking iterations (black markers): From left-to-right: a)
tracking efficiency, b) fake rate, and c) duplicate rate.

Figure 1 shows comparisons of basic physics performance markers. Tracking
efficiency is comparable overall; efficiency vs. $\eta$ (not shown) indicates
small gains in the endcap region ($2.4<|\eta|<2.8$). Fake rate is improved
overall with reduction improving with increasing $|\eta|$. Duplicate rate is
slightly increased but has been subsequently improved with further iteration-
specific tuning of the duplicate removal algorithm.

Figure 2: Comparison of computational performance of legacy track
reconstruction (red markers) and the new Run3 tracking configuration where
mkFit is used for 6 of twelve tracking iterations (black markers). Plots show
relative times of tracking steps for (left) iterations that use mkFit and
(right) all iterations.

Vectorization and threading scaling tests for initial iteration imply that,
according to Amdahlâs Law, ~70% of operations are vectorized and that more
than 95% of code is effectively parallelized. Computational speedups when
using mkFit are shown in figure 2. For all iterations where mkFit is used, the
observed track building time is reduced by ~3.5x (the best observed reduction
for one mkFit iteration is 6.7x). Note that track building with mkFit takes
less time than seeding, and about the same time as the final fit.

When all iterations, including non-mkFit based ones, are considered the
building time is reduced by ~1.7x. This translates to a 25% reduction of total
tracking time (including seeding, and final fit) and, overall, results in a
10-15% increase of Run3 reconstruction job event throughput.

##  3 Generalizations for iterative tracking & HL-LHC

CMSSW is a multi-threaded, module-based event processing framework that
instantiates and runs modules as well as manages data sources (both event data
and longer-lived data products) according to the dependencies generated by the
modules themselves when the job configuration is processed. As such, each
module can be instantiated multiple times and associated with different
configurations and data-sources, including parallel processing of several
events. This requires complete separability of configuration and module
instance state as well as complete absence of any non-constant global state.
Further, in the mkFit case, each iteration has its own set of parameters that
control and steer the functioning of core tracking algorithms as well as
separate implementations of standard functions, as mentioned in section 2.1.

While basic, algorithmic modifications had to be made to make mkFit conform to
the iterative tracking of CMSSW — i.e., to support forward search, backward
fit, and backward search — the majority of these changes amounted to
generalizations of the algorithms, along with mechanisms for expressing
different modes of behavior through configuration structures and intermediate-
level code that steers the algorithms. The previous section dealt with how the
code is structured and how it operates; this section addresses the design of
configuration structures and associated processing that allows the core of the
code, which is independent of experiment and geometry, to run in accordance
with the given detector description, algorithm tuning, and required standard
functions.

###  3.1 Geometry & detector description

In mkFit terminology a _layer_ denotes an $r-z$ bounding box in global
cylindrical space where hits belonging to the said layer are expected to be
found. It usually corresponds in some way to detector construction or readout
layers, but it does not have to:444E.g., mono and stereo hits from the same
silicon-strip detector layer in CMS are split into separate mkFit layers. its
main purposes are to aggregate the hits, provide an easy way to specify layer
crossing sequences for each tracking region (called a _layer plan_), and allow
track search to proceed uniformly among a set of tracks. This reduces
complexity and allows for the vectorization of certain key computations,
including track candidate propagation, hit selection, and Kalman filter
calculations.

Prior to CMSSW integration, logically dividing CMS into nested layers was
sufficient to allow mkFit to roughly reproduce the physics performance of
CMSSW legacy tracking. However, as mkFit was being considered a drop-in
replacement for the existing tracking implementation, additional detector-
module identification had to be included in mkFit’s geometry description to
enable it to pick up multiple hits from overlapping modules within the same
layer. Moreover, to support the Phase-2 upgrade geometry, which includes
axially tilted detector modules, further information had to be provided
(module position, normal and $\varphi$-direction vectors).

The layer boundaries, module details, and material properties are all
extracted during the CMSSW job setup by traversing all inner tracker modules.
For standalone usage this information gets exported into a binary file.

###  3.2 Configuration structures

As mentioned in the introduction to this section, mkFit code needs to run
concurrently within the main process, where each execution module is
configured for its specific tracking iteration. As there can be no static or
global data, the required configuration (or the relevant fragments) needs to
be passed down the execution stack or stored in local objects. It is therefore
important that the configuration data are structured in a way to facilitate
such usage.

The top-level configuration for each tracking iteration is represented by the
class IterationConfig. It contains flags that control which steps of the track
finding algorithm (see section 2.2) need to be performed, the standard
functions that are to be used for this iteration (described in more detail in
the next subsection), some high-level parameters for seed and duplicate
cleaning, and the following structures.

  * •

Layer traversal plans for all tracking regions.

  * •

Tracking parameters (e.g., maximum number of missed layers, $\chi^{2}$ cuts,
quality filter parameters) encapsulated in class IterationParams, with two
separate instances for forward and backward search.

  * •

Iteration-specific layer information, stored in class IterationLayerConfig,
which holds parameters guiding hit search and selection algorithms. These are
stored in a vector, with one instance per layer.

The CMSSW module system is typically configured via Python scripts; this
requires a rather tight coupling at the level of C++ code to parse incoming
data. As the above mkFit configuration is rather elaborate, it was accepted as
a compromise that all mkFit configuration can be loaded from (and saved into)
JSON files. Each iteration’s configuration is stored in a separate file
(stored as a part of CMSSW release) and the name of this file is then passed
to the mkFit CMSSW module during instantiation.

To allow for an easy modification of a small number of parameters, reading of
partial JSON overrides is fully supported: the default base configuration is
read from the CMSSW release and then existing in-memory representation gets
patched or overridden via simple additional JSON files or strings. Some
frequently used parameters can also be set via the Python interface, e.g., to
tune mkFit performance for heavy-ion operations.

Plugin-style configuration is still supported in standalone mode and is, in
fact, used to generate the default JSON files for the CMSSW operation.

###  3.3 Standard function catalogs

While adding support for multiple iterations and for Phase-2 tracking it
became obvious that using a single standard function and putting additional
parameters into IterationConfig structure does not scale and that a more
flexible configuration mechanism for standard functions is required for the
following tasks:

  * •

seed cleaning & partitioning â defined per iteration;

  * •

candidate filters, pre- and post-backward fit â defined per iteration;

  * •

duplicate cleaning â defined per iteration; and

  * •

candidate scoring â defined per iteration with a possible override for each
tracking region.

To provide a mechanism for registering different, specialized implementations
of these standard functions, and to be able to choose them at configuration
time, _standard function catalogs_ have been introduced. For each type of
function above, a thread-safe catalog with string keys and std::function value
type is provided. The catalogs are populated via static object initializers in
source files that contain the standard function codes. As std::function
objects are exported, the functions themselves can be hidden in anonymous
namespaces. Further, function templates can be used to inject compile-time
parameters and, in simple cases, the registered functions can be direct lambda
expressions.

With this infrastructure in place, JSON files can simply specify the names
(strings) associated in the catalog with the desired function. After
configuration loading and setup is complete the names get resolved into
std::function objects for fast access and become available through the
IterationConfig structure.

##  4 Ongoing & future work

At present, the described changes are being used to further tune Phase-1 CMS
iterations. Algorithmic improvements in the processing of layers and multi-
layer scoring are being investigated, with the goal of extending the usage of
mkFit beyond the current five tracking iterations, as well as improving
computational performance for currently supported use-cases. In parallel,
Phase-2 tracking is being developed, currently still focusing on a single,
initial tracking iteration, while the specifics of track propagation and
Kalman updates required for the support of tilted modules are being worked
out.

The final fit is now the most time-consuming tracking task in iterations using
mkFit. With the latest additions to geometry description, it should be
feasible to effectively use mkFit for this task as well, and we are
investigating the required developments in this area, along with possible
improvements to existing backward-fit and backward-search algorithms.

There has been a recent intense development of CMS Phase-2 Line Segment
Tracking (LST) lst-for-hllhc ; lst-cms-dpnote , a highly parallelizable
algorithm that can run efficiently on GPUs, which is showing great promise for
both offline and high-level trigger usage. We are planning to explore possible
synergetic development with the LST project, aiming for a hybrid approach
where LST performs the initial, fast track finding and mkFit provides final
steps such as backward fit, overlap hit search, and final fitting.

##  5 Conclusion

mkFit is in production mode for the CMS experiment since Run3 of the LHC, used
as a drop-in replacement for the legacy tracking code for five out of twelve
iterations, with equivalent physics performance and with overall tracking time
reduction of ~25%. Work has started to support CMS Phase-2 tracking geometry
and events with increased pileup, where some of required changes are described
in this paper, namely: generalizations of geometry description, multi-
iteration configuration, and introduction of catalogs of standard functions.
Exploration of extending mkFit to also cover the final-fit and to operate
synergistically with other track finding algorithms is in progress.

### Acknowledgements

This work was supported by the U.S. National Science Foundation under
Cooperative Agreements OAC-1836650 and PHY-2121686 and grant NSF-PHY-1912813.

## References

  * (1) R. Frühwirth, _Application of Kalman filtering to track and vertex fitting,_ Nucl. Instrum. Meth. A262, 440â450 (1987)
  * (2) CMS Collaboration, _The CMS experiment at the CERN LHC_ JINST 3, S08004 (2008)
  * (3) CMS Collaboration, _Development of the CMS detector for the CERN LHC Run 3,_ arXiv:2309.05466 [physics.ins-det] (2023) https://arxiv.org/abs/2309.05466
  * (4) S. Lantz _et al._ , _Speeding up particle track reconstruction using a parallel Kalman filter algorithm,_ JINST 15, P09030 (2020)
  * (5) G. Cerati _et al._ , _Parallelized and Vectorized Tracking Using Kalman Filters with CMS Detector Geometry and Events,_ EPJ Web of Conferences 214, 02002 (2019)
  * (6) CMS Collaboration, _Performance of Run 3 track reconstruction with the mkFit algorithm,_ CERN-CMS-DP-2022-018 (2022) https://cds.cern.ch/record/2814000
  * (7) CMS collaboration, _Description and performance of track and primary-vertex reconstruction with the CMS tracker,_ JINST 9, P10009 (2014) [arXiv:1405.6569]
  * (8) P. Chang _et al._ , _Line Segment Tracking in the High-luminosity LHC,_ in these proceedings
  * (9) CMS Collaboration, _Performance of Line Segment Tracking algorithm at HL-LHC,_ CERN-CMS-DP-2023-019 (2023) https://cds.cern.ch/record/2857438