$\displaystyle\|\widehat{\mathbf{W}}\mathbf{\Lambda}^{*\frac{1}{2}}\widehat{\mathbf{W}}^{\top}-\widehat{\mathbf{\Lambda}}^{\frac{1}{2}}\|_{2}\leq
C\frac{\lambda_{1}^{2}(\mathbf{\Sigma}^{*})}{\lambda_{r}^{3/2}(\mathbf{L}^{*})}\sqrt{\frac{\tau\vee\log
p}{n}}.$ (16)

Now we can plug the upper bounds for
$\|\widehat{\mathbf{U}}\widehat{\mathbf{W}}-\mathbf{U}^{*}\|_{2}$,
$\|\widehat{\mathbf{U}}\widehat{\mathbf{W}}-\mathbf{U}^{*}\|_{2\rightarrow\infty}$,
and
$\|\widehat{\mathbf{\Lambda}}^{\frac{1}{2}}-\widehat{\mathbf{W}}\mathbf{\Lambda}^{*\frac{1}{2}}\mathbf{W}^{\top}\|_{2}$
into (14) and the proof is complete. ∎

###### Proof of Lemma 2.

Define function
$h(x)=\sum_{i,j}\frac{u_{i}u_{j}}{v_{i}+v_{j}}x^{v_{i}+v_{j}}$, then some
calculation shows that

$\frac{dh(x)}{dx}=\sum_{i,j}u_{i}u_{j}x^{v_{i}+v_{j}-1}=\frac{1}{x}(\sum_{i}u_{i}x^{v_{i}})^{2},$

which implies $\frac{dh(x)}{dx}\geq 0$ for $x>0$. Since $h(x)$ is continuous
on $[0,\infty)$, and differentiable on $(0,\infty)$, Taylor’s theorem suggests
that $\sum_{i,j}\frac{u_{i}u_{j}}{v_{i}+v_{j}}=h(1)\geq h(0)=0$.

∎

##  Appendix E Additional Simulation Studies

In this section, we present additional simulation studies. As in the main
text, we evaluate the low-rank graph quilting approaches based on the
Frobenius norm between the imputed and true covariance matrices, as well as
F-1 scores from the resulting graphical Lasso estimates and selected edge set
from the MADgq approach with respect to the graphical Lasso estimate using the
full covariance matrix. In Tables 5(a) and 5(b), we show results for changing
the number of observations generated during the simulation procedure,
specifically with $n=1000,1500,2000,2500,$ and $3000$ for $p=100$ and with
$K=2$ patches of $o=60$ features. We see that the average Frobenius norm
between the imputed and true full covariance matrices decreases with a larger
number of observations. However, the F-1 scores for edge selection do not
fluctuate substantially with changing the number of observations, though the
standard deviation of the estimates decreases with increasing the number
samples. We also show the performance of the low-rank graph quilting methods
when the low-rankness of the simulated full covariance matrix changes. We
measure this by the ratio of the first eigenvalue of the covariance matrix to
the last; specifically, for $n=2000$ and $p=100$, with $K=2$ patches of size
$o=60$, we look at $\lambda_{1}/\lambda_{p}=500,200,50,10,$ and $3$; here, a
smaller ratio represents a weaker low-rank structure. The Frobenius norm of
the imputed covariance matrices decreases when the leading eigenvalue is
smaller, though this may be an impact of the generated full covariance
matrices rather than the estimation procedure. Also, as we would expect, when
the true covariance matrix has less of a clear rank 1 structure, all of the
low-rank graph quilting methods perform worse in terms of edge selection for
the resulting graphical model estimate. On the other hand, the MADgq method
performs much better. This matches with the assumptions made by the two
different types of methods.

Model | $\bm{n=1000}$ | $\bm{n=1500}$ | $\bm{n=2000}$ | $\bm{n=2500}$ | $\bm{n=3000}$
---|---|---|---|---|---
BSVDgq, LR | 0.836 (0.091) | 0.852 (0.068) | 0.716 (0.157) | 0.678 (0.124) | 0.578 (0.090)
BSVDgq, ALR | 0.761 (0.060) | 0.790 (0.144) | 0.735 (0.122) | 0.724 (0.129) | 0.645 (0.150)
LRFgq, LR | 0.920 (0.210) | 0.968 (0.195) | 0.844 (0.186) | 0.834 (0.159) | 0.843 (0.202)
LRFgq, ALR | 1.015 (0.339) | 0.759 (0.140) | 0.715 (0.142) | 0.786 (0.096) | 0.667 (0.201)
NNgq, LR | 0.831 (0.140) | 0.787 (0.130) | 0.699 (0.158) | 0.661 (0.143) | 0.576 (0.135)
NNgq, ALR |  0.503 (0.094) |  0.465 (0.093) |  0.408 (0.097) |  0.371 (0.065) |  0.363 (0.058)
Zero Imputation | 5.679 (0.530) | 5.736 (0.592) | 5.976 (0.597) | 5.790 (0.458) | 5.592 (0.650)
Model | $\bm{\lambda_{1}/\lambda_{p}\approx 500}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 200}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 50}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 10}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 3}$
BSVDgq, LR | 0.716 (0.157) | 0.657 (0.095) | 0.516 (0.069) | 0.486 (0.058) | 0.311 (0.036)
BSVDgq, ALR | 0.735 (0.122) | 0.636 (0.095) | 0.514 (0.091) | 0.470 (0.043) | 0.308 (0.024)
LRFgq, LR | 0.844 (0.186) | 0.659 (0.099) | 0.529 (0.052) | 0.482 (0.026) | 0.317 (0.035)
LRFgq, ALR | 0.715 (0.142) | 0.646 (0.124) | 0.545 (0.065) | 0.454 (0.057) | 0.284 (0.025)
NNgq, LR | 0.699 (0.158) | 0.409 (0.062) | 0.380 (0.035) | 0.283 (0.034) | 0.068 (0.022)
NNgq, ALR |  0.408 (0.097) |  0.172 (0.037) |  0.184 (0.046) |  0.149 (0.036) |  0.064 (0.026)
Zero Imputation | 5.976 (0.597) | 1.991 (0.234) | 0.603 (0.140) | 0.283 (0.035) | 0.069 (0.024)

(a) Frobenius norms for the imputed covariance matrices from the low-rank
graph quilting, the MADgq procedure, and zero imputation methods as compared
to the full covariance matrix from simulated data. Best performing are
boldfaced.

Model | $\bm{n=1000}$ | $\bm{n=1500}$ | $\bm{n=2000}$ | $\bm{n=2500}$ | $\bm{n=3000}$
---|---|---|---|---|---
BSVDgq, LR | 0.845 (0.049) | 0.840 (0.048) | 0.856 (0.045) | 0.862 (0.044) |  0.858 (0.035)
BSVDgq, ALR | 0.843 (0.037) | 0.795 (0.041) | 0.855 (0.039) | 0.861 (0.038) | 0.857 (0.029)
LRFgq, LR | 0.844 (0.044) |  0.841 (0.037) |  0.856 (0.041) | 0.863 (0.045) | 0.833 (0.039)
LRFgq, ALR |  0.867 (0.031) | 0.830 (0.038) | 0.854 (0.040) |  0.866 (0.030) | 0.833 (0.041)
NNgq, LR | 0.859 (0.041) | 0.855 (0.039) | 0.846 (0.037) | 0.845 (0.034) | 0.828 (0.038)
NNgq, ALR | 0.462 (0.054) | 0.533 (0.042) | 0.739 (0.039) | 0.731 (0.042) | 0.738 (0.041)
MADgq | 0.287 (0.056) | 0.343 (0.042) | 0.314 (0.051) | 0.332 (0.045) | 0.344 (0.042)
Zero Imputation | 0.333 (0.042) | 0.323 (0.041) | 0.314 (0.054) | 0.309 (0.043) | 0.315 (0.049)
Model | $\bm{\lambda_{1}/\lambda_{p}\approx 500}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 200}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 50}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 10}$ | $\bm{\lambda_{1}/\lambda_{p}\approx 3}$
BSVDgq, LR | 0.856 (0.045) | 0.849 (0.044) | 0.828 (0.042) | 0.615 (0.049) | 0.571 (0.055)
BSVDgq, ALR | 0.854 (0.039) | 0.754 (0.041) | 0.773 (0.037) |  0.718 (0.044) | 0.574 (0.048)
LRFgq, LR |  0.856 (0.041) |  0.848 (0.038) |  0.829 (0.044) | 0.614 (0.038) | 0.581 (0.053)
LRFgq, ALR | 0.855 (0.040) | 0.821 (0.035) | 0.773 (0.041) | 0.662 (0.039) | 0.584 (0.052)
NNgq, LR | 0.846 (0.037) | 0.838 (0.031) | 0.822 (0.044) | 0.465 (0.028) | 0.541 (0.050)
NNgq, ALR | 0.739 (0.039) | 0.752 (0.041) | 0.752 (0.039) | 0.413 (0.035) | 0.286 (0.044)
MADgq | 0.317 (0.031) | 0.322 (0.039) | 0.363 (0.044) | 0.517 (0.045) |  0.689 (0.049)
Zero Imputation | 0.314 (0.033) | 0.307 (0.034) | 0.345 (0.032) | 0.468 (0.034) | 0.424 (0.036)

(b) F-1 scores for edge selection from the low-rank graph quilting methods,
the MADgq procedure, and zero imputation, when compared the edges selected
from the graphical Lasso on the full empirical covariance matrix.

Table 5: Results from additional simulation studies. Best performing methods
are boldfaced.

##  Appendix F Neuronal Tuning Results

In Table 6, we compare the selected functional connections by each low-rank
graph quilting method in the unobserved portion of the covariance matrix by
neuronal tuning. As mentioned in the main text, we expect Neurons are
categorized by 8 categories of angular tuning and 5 categories of frequency
tuning. Generally, the performance of the different low-rank graph quilting
methods is relatively similar for finding functional connections with matching
angular tuning, while the approximate low-rank assumption appears to generally
perform better for frequency tuning. We also see that all of the low-rank
graph quilting methods find functional connections between neurons with the
same tuning at similar rates as when the full data is observed and at a
substantially higher rate compared the MADgq method.

Model | Angular Tuning | Frequency Tuning
---|---|---
(Full Data) | 0.208 | 0.357
BSVDgq, LR | 0.247 (0.057) | 0.353 (0.039)
BSVDgq, ALR | 0.241 (0.054) | 0.413 (0.048)
LRFgq, LR | 0.226 (0.044) | 0.319 (0.045)
LRFgq, ALR | 0.251 (0.056) | 0.392 (0.051)
NNgq, LR | 0.284 (0.063) | 0.452 (0.061)
NNgq, ALR | 0.257 (0.052) | 0.377 (0.042)
MADgq | 0.145 (0.027) | 0.251 (0.040)

Table 6: Proportion of edges from the imputed portion of the estimated
functional connectivity network which link two neurons of the same tuning
category for each low-rank graph quilting method. Results for $K=2,o=150$ in
the Allen Institute data.

##  Appendix G Real Data Inspired Simulation

(a) Allen Institute data.

(b) Simulation data.

(c) Imputed NNgq.

(d) Imputed BSVDgq.

Figure 4: Eigendecompositions on empirical and imputed covariance matrices.

As mentioned in the main text, we see from our empirical studies that the
gradient descent and block singular value decomposition methods tend to
produce higher F-1 scores in the simulated data, while the nuclear norm
penalization methods are more likely to have higher F-1 scores for both of the
real data experiments. This can be explained by the structure of the true
underlying covariance matrices as well as those of the imputed covariance
matrices from each method. Figures 4(a) and 4(b) show the eigenvalues from the
decomposition on the covariance matrix for the Allen Institute data and from
the simulated data. As we can see, the real data have exponentially decaying
eigenvalues, while the simulated data look almost exactly rank 1, with all
other eigenvalues near 0. In Figures 4(c) and 4(d), we see eigenvalues from
the decomposition of the imputed covariance matrices from the nuclear norm
penalization and the block singular value decomposition methods. In the former
we see an exponentially decaying eigenvalue sequence, and in the latter, we
see more of an exact low-rank structure. We explore this further in the
simulation study below.

Model | $\bm{n=10000}$ | $\bm{n=15000}$ | $\bm{n=20000}$ | $\bm{n=25000}$ | $\bm{n=30000}$
---|---|---|---|---|---
BSVDgq, LR | 3.258 (0.431) | 3.115 (0.450) | 3.215 (0.414) | 2.984 (0.423) | 3.152 (0.397)
BSVDgq, ALR | 3.249 (0.459) | 3.076 (0.417) | 3.212 (0.458) | 2.976 (0.426) | 3.226 (0.384)
LRFgq, LR | 3.365 (0.446) | 3.080 (0.433) | 3.181 (0.439) | 3.035 (0.471) | 3.167 (0.404)
LRFgq, ALR | 3.376 (0.395) | 3.092 (0.428) | 3.220 (0.465) | 3.055 (0.467) | 3.098 (0.435)
NNgq, LR | 2.697 (0.340) | 2.705 (0.341) | 2.578 (0.397) | 2.941 (0.387) | 2.768 (0.335)
NNgq, ALR |  1.405 (0.301) |  1.634 (0.340) |  1.418 (0.317) |  1.105 (0.305) |  1.302 (0.265)
Zero Imputation | 4.213 (0.515) | 4.097 (0.490) | 4.063 (0.475) | 4.244 (0.506) | 4.221 (0.464)

(a) Frobenius norms for the imputed covariance matrices from the low-rank
graph quilting and zero imputation methods as compared to the full covariance
matrix.

Model | $\bm{n=10000}$ | $\bm{n=15000}$ | $\bm{n=20000}$ | $\bm{n=25000}$ | $\bm{n=30000}$
---|---|---|---|---|---
BSVDgq, LR | 0.795 (0.042) | 0.793 (0.037) | 0.798 (0.035) | 0.791 (0.033) | 0.796 (0.034)
BSVDgq, ALR | 0.866 (0.037) | 0.810 (0.036) | 0.825 (0.029) | 0.838 (0.030) | 0.832 (0.029)
LRFgq, LR | 0.791 (0.039) | 0.794 (0.035) | 0.787 (0.033) | 0.789 (0.032) | 0.793 (0.028)
LRFgq, ALR | 0.860 (0.034) | 0.807 (0.035) | 0.824 (0.031) | 0.838 (0.027) | 0.832 (0.026)
NNgq, LR | 0.805 (0.033) | 0.764 (0.038) | 0.788 (0.030) | 0.715 (0.027) | 0.711 (0.028)
NNgq, ALR |  0.875 (0.031) |  0.866 (0.026) |  0.887 (0.023) |  0.853 (0.023) |  0.860 (0.024)
MADgq | 0.442 (0.050) | 0.449 (0.045) | 0.401 (0.049) | 0.461 (0.044) | 0.445 (0.042)
Zero Imputation | 0.558 (0.048) | 0.411 (0.049) | 0.541 (0.046) | 0.559 (0.042) | 0.508 (0.044)

(b) F-1 scores for edge selection for the low-rank graph quilting methods
compared the edges selected from the graphical Lasso on the fully observed
covariance.

Table 7: Results from the low-rank graph quilting methods, the MADgq
procedure, and zero imputation compared to results for the full empirical
covariance, in new data simulated using the empirical covariance matrix from
the Allen Institute data. The sample size $n$ is changed.

Here, instead of creating a rank-1 covariance matrix as the true fully
observed covariance matrix in Section 3, we instead use the empirical
covariance matrix from the Allen Institute calcium imaging data from Section
4. Specifically, we take the covariance matrix from the real data as the true
underlying covariance matrix, and use it to simulate new data from a
multivariate Gaussian distribution. As done with the previous simulations in
Section 3, we then create from this simulated data a partially observed
covariance matrix with $k$ patches of $o$ features each. With these, we impute
the covariance matrices using the low-rank graph quilting methods and zero
imputation, then apply the graphical Lasso to them to obtain graphical model
estimates. Additionally, we apply the MADgq to the partially observed
covariance matrix to estimate the edges of the graphical model directly. The
Frobenius norm and F-1 score metrics are used to compare the methods. As
above, hyperparameter values for the nuclear norm and low-rank factorization
covariance imputation methods are selected by optimal Frobenius norm scores
with respect to the true full covariance matrix, while hyperparameter
selection for the graphical Lasso step, when applicable, is performed using
oracle sparsity tuning with respect to the full empirical covariance matrix.

Model | $\bm{K=2,o=130}$ | $\bm{K=2,o=140}$ | $\bm{K=2,o=150}$ | $\bm{K=2,o=160}$ | $\bm{K=2,o=170}$
---|---|---|---|---|---
BSVDgq, LR | 4.014 (0.413) | 3.552 (0.408) | 3.215 (0.414) | 2.228 (0.397) | 2.277 (0.347)
BSVDgq, ALR |  3.617 (0.469) | 3.521 (0.404) | 3.212 (0.458) | 2.267 (0.360) | 2.255 (0.378)
LRFgq, LR | 3.966 (0.449) | 3.501 (0.467) | 3.181 (0.439) | 2.220 (0.356) | 2.075 (0.305)
LRFgq, ALR | 3.628 (0.454) | 3.503 (0.411) | 3.220 (0.465) | 2.261 (0.394) | 2.281 (0.354)
NNgq, LR | 5.137 (0.445) | 4.113 (0.385) | 2.578 (0.397) | 1.395 (0.166) | 1.282 (0.217)
NNgq, ALR | 4.489 (0.401) |  1.677 (0.316) |  1.418 (0.317) |  0.807 (0.087) |  0.821 (0.095)
Zero Imputation | 5.138 (0.467) | 4.563 (0.447) | 4.063 (0.475) | 3.149 (0.409) | 3.073 (0.423)

(a) Frobenius norms for the imputed covariance matrices from the low-rank
graph quilting and zero imputation methods as compared to the full covariance
matrix.

Model | $\bm{K=2,o=130}$ | $\bm{K=2,o=140}$ | $\bm{K=2,o=150}$ | $\bm{K=2,o=160}$ | $\bm{K=2,o=170}$
---|---|---|---|---|---
BSVDgq, LR | 0.612 (0.040) | 0.693 (0.036) | 0.798 (0.035) | 0.820 (0.035) | 0.845 (0.026)
BSVDgq, ALR | 0.660 (0.038) | 0.765 (0.033) | 0.825 (0.029) | 0.871 (0.024) | 0.903 (0.023)
LRFgq, LR | 0.608 (0.041) | 0.701 (0.033) | 0.787 (0.033) | 0.818 (0.027) | 0.843 (0.033)
LRFgq, ALR |  0.661 (0.038) | 0.766 (0.031) | 0.824 (0.036) | 0.865 (0.029) | 0.897 (0.021)
NNgq, LR | 0.422 (0.038) | 0.525 (0.034) | 0.788 (0.030) | 0.824 (0.030) | 0.819 (0.029)
NNgq, ALR | 0.550 (0.039) |  0.783 (0.031) |  0.887 (0.023) |  0.884 (0.026) |  0.905 (0.025)
MADgq | 0.336 (0.050) | 0.365 (0.046) | 0.401 (0.049) | 0.436 (0.045) | 0.592 (0.044)
Zero Imputation | 0.423 (0.048) | 0.525 (0.049) | 0.541 (0.046) | 0.515 (0.050) | 0.596 (0.041)

(b) F-1 scores for edge selection for the low-rank graph quilting methods
compared the edges selected from the graphical Lasso on the fully observed
covariance.

Table 8: Results from the low-rank graph quilting methods, the MADgq
procedure, and zero imputation compared to results for the full empirical
covariance, in new data simulated using the empirical covariance matrix from
the Allen Institute data. The patch size $o$ is changed.

Model | $\bm{K=2,o=150}$ | $\bm{K=3,o=100}$ | $\bm{K=4,o=75}$ | $\bm{K=5,o=60}$ | $\bm{K=6,o=50}$
---|---|---|---|---|---
BSVDgq, LR | 3.215 (0.414) | 3.513 (0.435) | 3.601 (0.464) | 5.302 (0.463) | 5.586 (0.524)
BSVDgq, ALR | 3.212 (0.458) | 3.574 (0.406) | 3.556 (0.475) | 3.628 (0.424) | 4.200 (0.506)
LRFgq, LR | 3.181 (0.439) | 3.923 (0.457) | 3.852 (0.425) | 5.010 (0.503) | 5.254 (0.467)
LRFgq, ALR | 3.220 (0.465) | 3.576 (0.445) | 3.557 (0.433) | 3.529 (0.479) | 4.202 (0.498)
NNgq, LR | 2.578 (0.397) | 4.138 (0.411) | 4.838 (0.424) | 4.849 (0.453) | 5.647 (0.459)
NNgq, ALR |  1.418 (0.317) |  2.192 (0.350) |  2.125 (0.363) |  2.986 (0.327) |  3.957 (0.386)
Zero Imputation | 4.063 (0.475) | 5.442 (0.457) | 5.282 (0.456) | 5.215 (0.440) | 5.734 (0.482)

(a) Frobenius norms for the imputed covariance matrices from the low-rank
graph quilting and zero imputation methods as compared to the full covariance
matrix.

Model | $\bm{K=2,o=150}$ | $\bm{K=3,o=100}$ | $\bm{K=4,o=75}$ | $\bm{K=5,o=60}$ | $\bm{K=6,o=50}$
---|---|---|---|---|---
BSVDgq, LR | 0.798 (0.035) | 0.736 (0.039) | 0.628 (0.032) | 0.608 (0.033) | 0.452 (0.045)
BSVDgq, ALR | 0.825 (0.029) | 0.757 (0.036) | 0.689 (0.043) | 0.660 (0.036) | 0.482 (0.049)
LRFgq, LR | 0.787 (0.033) | 0.739 (0.038) | 0.649 (0.045) | 0.610 (0.043) | 0.463 (0.042)
LRFgq, ALR | 0.824 (0.031) | 0.756 (0.032) | 0.682 (0.038) |  0.673 (0.035) |  0.499 (0.048)
NNgq, LR | 0.788 (0.030) | 0.563 (0.040) | 0.578 (0.042) | 0.535 (0.039) | 0.353 (0.044)
NNgq, ALR |  0.887 (0.023) |  0.789 (0.032) |  0.690 (0.046) | 0.545 (0.042) | 0.412 (0.042)
MADgq | 0.401 (0.049) | 0.363 (0.045) | 0.336 (0.048) | 0.242 (0.039) | 0.197 (0.031)
Zero Imputation | 0.541 (0.046) | 0.456 (0.044) | 0.458 (0.039) | 0.443 (0.043) | 0.257 (0.028)

(b) F-1 scores for edge selection for the low-rank graph quilting methods
compared the edges selected from the graphical Lasso on the fully observed
covariance.

Table 9: Results from the low-rank graph quilting methods, the MADgq
procedure, and zero imputation compared to results for the full empirical
covariance, in new data simulated using the empirical covariance matrix from
the Allen Institute data. The number of patches $k$ is changed.

In this simulation, we let $p=227$ as in the Allen Institute dataset. Tables
7(a) and 7(b) show results for $n=10000,15000,20000,25000,$ and $30000$, with
$k=2$ and $o=150$. Tables 8(a) and 8(b) for $o=130,140,150,160$, and $170$
with $k=2$ and $n=20000$. Tables 9(a) and 9(b) show the results for
$k=2,3,4,5$ and $6$, for $n=20000$ and $o=150,100,75,60$, and $50$,
respectively. In these simulations, we see that the nuclear norm penalization
methods now generally have higher F-1 scores compared the block singular value
decomposition and the gradient descent methods for edge recovery of the graph
from using the full covariance matrix; this falls more in line with what was
seen in the results from the real data studies in Section 4, and the opposite
of the results seen in Section 3. Otherwise, we see similar patterns in terms
of changes in performance based on changes to simulation parameters that we do
in Section 3.


# Multilevel Method for Thermal Radiative Transfer Problems with Method of
Long Characteristics for the Boltzmann Transport Equation

Joseph M. Coale<EMAIL_ADDRESS>Dmitriy Y. Anistratov<EMAIL_ADDRESS>Los Alamos National Laboratory, Los Alamos, NM 87545  Department of Nuclear
Engineering, North Carolina State University, Raleigh, NC 27695

###### Abstract

In this paper analysis is performed on a computational method for thermal
radiative transfer (TRT) problems based on the multilevel quasidiffusion
(variable Eddington factor) method with the method of long characteristics
(ray tracing) for the Boltzmann transport equation (BTE). The method is
formulated with a multilevel set of moment equations of the BTE which are
coupled to the material energy balance (MEB). The moment equations are exactly
closed via the Eddington tensor defined by the BTE solution. Two discrete
spatial meshes are defined: a material grid on which the MEB and low-order
moment equations are discretized, and a grid of characteristics for solving
the BTE. Numerical testing of the method is completed on the well-known Fleck-
Cummings test problem which models a supersonic radiation wave propagation.
Mesh refinement studies are performed on each of the two spatial grids
independently, holding one mesh width constant while refining the other. We
also present the data on convergence of iterations.

###### keywords:

Boltzmann transport equation, high-energy density physics, thermal radiative
transfer long characteristics, ray tracing, multilevel methods, quasidiffusion
method, variable Eddington factor.

##  1 Introduction

In this paper we investigate the performance of computational methods for
high-energy density physics (HEDP) simulations [1, 2]. Phenomena that occur in
the high-energy density (HED) regime are typically modeled with complex
multiphysical systems of partial differential equations (PDEs). Energy
redistribution in the HED regime is dominated by the mechanism of radiative
transfer, where energy is transported by the emission and absorption of photon
radiation. The radiative transfer equation is the Boltzmann transport equation
(BTE) which models how photons propagate through and interact with matter.
Thus, the BTE is an essential component in the models of HED phenomena. These
systems of PDEs exhibit several properties which present challenges for their
numerical simulation, including: (i) strong nonlinearity, (ii) tight coupling
between equations, (iii) multiscale behavior in space-time, (iv) high-
dimensionality.

Discretization of the equations involved in HEDP simulations must be handled
with these properties in mind to preserve the fundamental physical behavior of
the modeled phenomena. To this end the method of long characteristics (MOLC),
a.k.a. ray-tracing schemes (RTS), for the BTE hold some advantage compared to
other discretization schemes. Numerical techniques based on MOLC/RTS for
solving particle transport problems are formulated by means of the integral
form of the BTE. They belong to the family of the method of characteristics
that has a rich and long history of development [3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16].

This terminology reflects that the MOLC solves the BTE along characteristics
extending over the entire spatial domain. The MOLC doesn’t employ
interpolation of solution on faces of cells. Since the solution is not
interpolated on the spatial mesh, discontinuities can be well resolved and
carried across the problem domain. The BTE solution will also be defined on a
unique discrete grid separate from the multiphysics equations it becomes
coupled to. It is therefore straightforward to refine or coarsen the discrete
BTE grid independently of the multiphysics equations, allowing for
computational resources to be more effectively managed.

We consider multi-dimensional, time dependent problems that capture the
essential challenges in HEDP simulation. Specifically, the computational
method is developed on and analyzed with the fundamental thermal radiative
transfer (TRT) problem, formulated by the multigroup BTE for photons

$\frac{1}{c}\frac{\partial I_{g}}{\partial
t}+\boldsymbol{\Omega}\cdot\boldsymbol{\nabla}I_{g}+\varkappa_{g}(T)I_{g}=\varkappa_{g}(T)B_{g}(T)\,,$
(1) $I_{g}|_{\boldsymbol{r}\in\partial\Gamma}=I_{g}^{\text{in}}\ \ \text{for}\
\ \boldsymbol{\Omega}\cdot\boldsymbol{n}_{\Gamma}<0,\quad
I_{g}|_{t=0}=I_{g}^{0}\,,$ (2) $\boldsymbol{r}\in\Gamma,\ t>0,\
\boldsymbol{\Omega}\in\mathcal{S},\ g=1,\dots,G$

coupled to the material energy balance (MEB) equation that describes
radiation-matter energy exchange

$\frac{\partial\varepsilon(T)}{\partial
t}={\sum_{g=1}^{G}\bigg{(}\int_{4\pi}I_{g}d\Omega-
B_{g}(T)\bigg{)}\varkappa_{g}(T)},\quad T|_{t=0}=T^{0}\,,$ (3)

where $\boldsymbol{r}$ is spatial position, $\boldsymbol{\Omega}$ is the
direction of particle motion, $g$ is the photon frequency group index, $c$ is
the speed of light, $\Gamma$ is the spatial domain, $\partial\Gamma$ is the
domain boundary, $\boldsymbol{n}_{\Gamma}$ is the unit normal to
$\partial\Gamma$ and $\mathcal{S}$ is the unit sphere.
$I_{g}(\boldsymbol{r},\boldsymbol{\Omega},t)$ is the group specific intensity
of radiation, $T(\boldsymbol{r},t)$ is the material temperature,
$\varkappa_{g}(\boldsymbol{r},t;T)$ is the group photon opacity,
$\varepsilon(\boldsymbol{r},t;T)$ is the material energy density and
$B_{g}(\boldsymbol{r},t;T)$ is the group Planck black-body distribution
function. This TRT problem captures all fundamental challenges and features
associated with multiphysical HEDP problems, and serves as a useful platform
for the analysis of computational methods aimed at HEDP simulation.

To solve the TRT problem we employ the multilevel quasidiffusion (QD), a.k.a.
Variable Eddington Factor (VEF), methodology [17, 18, 19, 20, 21, 22, 23, 24,
25], which is essentially a nonlinear method of moments. The multilevel QD
(MLQD) method is defined by a system consisting of (i) the high-order BTE and
(ii) a hierarchy of low-order equations for moments of the radiation
intensity. The low-order equations are exactly closed through the Eddington
tensor and other linear-fractional factors that are weakly dependent on the
BTE solution. Multiphysical equations (i.e. the MEB) are coupled to the low-
order equations on the scale of the multiphysics, effectively reducing the
dimensionality of the problem. The MOLC/RTS is used to compute the high-order
solution of the BTE that is used to compute the Eddington tensor to close the
moment equations. The important property of the MLQD method is that the
hierarchy of high-order BTE and moment equations can be discretized
independently [26]. This enables one to take advantage of essentially
different type of discretization for high-order and low-order equations to
find an improvement in the numerical solution and computational efficiency.

This paper aims to address several fundamental questions regarding the use of
MOLC/RTS as part of the MLQD hierarchy of equations; Specifically as
generating the Eddington tensor closures for HEDP applications. The
discretization of the system of high-order and low-order equations are
performed on two grids. The low-order equations for moments of the intensity
and the MEB equation for temperature are approximated on the underlying mesh
defined by material properties. We refer to it as the spatial material mesh.
The ray-tracing technique in MOLC/RTS defines a mesh of characteristics over
the whole spatial domain. The characteristics stretch from one domain boundary
to another. This procedure also creates a subgrid of characteristics in each
spatial material cell. Note that the mesh of characteristics and the resulting
subgrid in spatial material cells depend on discrete directions. It can also
vary for different frequency groups. We investigate the essential relation
between the two distinct discrete grids by performing studies on the effects
of refining each grid independently of one another. This process is meant to
shed light on the relative importance of each individual grid and how much
effort should be given to refine either the characteristic mesh or underlying
material mesh to optimize solution accuracy and performance. Furthermore, the
impact of mesh sizing on iterative convergence of the solution is examined.

The remainder of the paper is organized as follows. In Sec. 2, the MLQD method
is formulated. In Sec. 3, we present details of a variant of MOLC/RTS for
solving the high-order BTE. The numerical results are presented in Sec. 4. We
conclude with a discussion in Sec. 5.

##  2 Multilevel QD/VEF Method for TRT

Two systems of moment equations of the BTE are constructed, the first of which
being the multigroup low-order quasidiffusion (LOQD) equations

$\frac{\partial E_{g}}{\partial
t}+\boldsymbol{\nabla}\cdot\boldsymbol{F}_{g}+c\varkappa_{g}(T)E_{g}=4\pi\varkappa_{g}(T)B_{g}(T)\,,$
(4a) $\frac{1}{c}\frac{\partial\boldsymbol{F}_{g}}{\partial
t}+c\boldsymbol{\nabla}\cdot(\boldsymbol{\mathfrak{f}}_{g}E_{g})+\varkappa_{g}(T)\boldsymbol{F}_{g}=0\,,$
(4b)

which solve for the multigroup radiation energy density
$E_{g}=\frac{1}{c}\int_{4\pi}I_{g}d\Omega$ and flux
$\boldsymbol{F}_{g}=\int_{4\pi}\boldsymbol{\Omega}I_{g}d\Omega$. The second
system of moment equations is the effective grey LOQD equations

$\frac{\partial E}{\partial
t}+\boldsymbol{\nabla}\cdot\boldsymbol{F}+c\langle\varkappa\rangle_{E}E=c\langle\varkappa\rangle_{B}a_{R}T^{4}\,,$
(5a) $\frac{1}{c}\frac{\partial\boldsymbol{F}}{\partial
t}+c\boldsymbol{\nabla}\cdot(\langle\boldsymbol{\mathfrak{f}}\rangle_{E}E)+\bar{\mathbf{K}}_{R}\boldsymbol{F}+\bar{\boldsymbol{\eta}}E=0\,,$
(5b)

which solve for the total radiation energy density $E=\sum_{g=1}^{G}E_{g}$ and
flux $\boldsymbol{F}=\sum_{g=1}^{G}\boldsymbol{F}_{g}$. The MEB equation is
cast in effective grey form

$\frac{\partial\varepsilon(T)}{\partial
t}=c\langle\varkappa\rangle_{E}E-c\langle\varkappa\rangle_{B}a_{R}T^{4},$ (6)

and coupled with Eqs. (5). The multigroup LOQD equations (4) are exactly
closed via the Eddington tensor computed by the solution of the BTE (1):

$\boldsymbol{\mathfrak{f}}_{g}=\frac{\int_{4\pi}\boldsymbol{\Omega}\otimes\boldsymbol{\Omega}{I}_{g}\
d\Omega}{\int_{4\pi}I_{g}d\Omega}\,.$ (7)

The spectrum average opacities and coefficients calculated by the solution of
multigroup LOQD equations (4) which define exact closures for the LOQD system
are

$\langle
u\rangle_{E}=\frac{\sum_{g=1}^{G}u_{g}E_{g}}{\sum_{g=1}^{G}E_{g}},\quad\bar{\mathbf{K}}_{R}=\text{diag}\big{(}\langle
u\rangle_{|F_{x}|},\langle u\rangle_{|F_{y}|},\langle
u\rangle_{|F_{z}|}\big{)},,$ (8) $\langle
u\rangle_{|F_{\alpha}|}=\frac{\sum_{g=1}^{G}u_{g}|F_{\alpha,g}|}{\sum_{g=1}^{G}|F_{\alpha,g}|},\quad\bar{\boldsymbol{\eta}}=\frac{\sum_{g=1}^{G}(\varkappa_{g}-\bar{\mathbf{K}}_{R})\boldsymbol{F}_{g}}{\sum_{g=1}^{G}E_{g}}\,.$
(9)

In sum the MLQD method for the TRT problem is formulated with:

  1. 1.

The BTE (1) discretized with the MOLC/RTS,

  2. 2.

The multigroup LOQD equations (Eqs. (4) & (7)) discretized with the second-
order finite volume (FV) scheme [27],

  3. 3.

The effective grey problem formed by the closed system of effective grey LOQD
equations and the MEB equation (Eqs. (5), (6), (8) & (9)). The spatial
approximation FV scheme for effective grey LOQD equations are algebraically
consistent with the scheme for the multigroup LOQD equations [27].

##  3 Method of Long Characteristics / Ray-Tracing Scheme

The MOLC/RTS is derived by first performing a change of coordinates from
$(x,y)$ to $(u,v)$ where $\boldsymbol{e}_{u}=\boldsymbol{\Omega}$ and
$\boldsymbol{e}_{u}\cdot\boldsymbol{e}_{v}=0$. The BTE discretized in time
with the backward-Euler scheme, along characteristics is given

$\frac{\partial I^{n}_{g}(u)}{\partial
u}+\tilde{\varkappa}^{n}_{g}(u)I_{g}^{n}(u)={Q}^{n}_{g}(u)\,,$ (10)
$\tilde{\varkappa}^{n}_{g}=\varkappa^{n}_{g}+\frac{1}{c\Delta t^{n}},\quad
Q^{n}_{g}=\varkappa_{g}^{n}B_{g}^{n}+\frac{1}{c\Delta t^{n}}I_{g}^{n-1}\,,$
(11)

and $\Delta t^{n}$ is the $n^{\text{th}}$ time step. Let all low-order
equations be discretized on an ‘underlying’ orthogonal material spatial grid.
The material temperature is a piece-wise function on the set material spatial
cells. As a result, the opacities and Planckian emission source are constant
in each material cell. With the MOLC/RTS, the BTE is discretized on a
characteristic grid traced over the underlying grid for $M$ discrete
directions $\boldsymbol{\Omega}_{m}$. For each discrete direction, there will
be $K_{m}$ characteristics traced over the spatial domain for a total of $K$
characteristics which formulate the high-order discrete grid. Let $u^{+}$ and
$u^{-}$ be the intersection (entry and exit) points of a given characteristic
with the boundaries of a single cell in the underlying grid. This forms a
segment of the characteristic with length $\ell=(u^{+}-u^{-})\sin(\zeta_{m})$
where $\zeta_{m}$ is the angle between $\boldsymbol{\Omega}_{m}$ and the $z$
axis. Integrating Eq. (10) along each characteristic segment gives

$\displaystyle
I_{k,s+1}=I_{k,s}e^{-\frac{\tilde{\varkappa}_{i}\ell_{k,s}}{\sin(\zeta)}}+\frac{{Q}_{i}}{\tilde{\varkappa}_{i}}\big{(}1-e^{-\frac{\tilde{\varkappa}_{i}\ell_{k,s}}{\sin(\zeta)}}\big{)},\quad{\bar{I}}_{k,s}=\alpha_{k,s}I_{k,s}+(1-\alpha_{k,s})I_{k,s+1}\,,$
(12) $\displaystyle k=1,\dots,K,\quad s=1,\dots,S_{k}\,,$

where $m$ and $g$ subscripts have been emitted for brevity. $I_{k,s+1}$ is the
radiation intensity at the outgoing face of the $s^{\text{th}}$ segment of the
$k^{\text{th}}$ characteristic, $\tilde{\varkappa}_{i}$ and ${Q}_{i}$ are the
modified opacity and source in the $i^{\text{th}}$ material cell that is
traced over by the $s^{\text{th}}$ segment. ${\bar{I}}_{k,s}$ is the segment-
average radiation intensity, the optical thickness of each segment is

$\tau_{k,s}=\frac{\tilde{\varkappa}_{i}\ell_{k,s}}{\sin(\zeta)},\quad\alpha_{k,s}=\frac{1}{\tau_{k,s}}-\frac{e^{-\tau_{k,s}}}{1-e^{-\tau_{k,s}}}\,.$
(13)

The area of a characteristic segment is

$A_{k,s}=\ell_{k,s}w_{k}\,,$ (14)

where $w_{k}$ is the width of the $k^{\text{th}}$ characteristic. Radiation
intensities on each characteristic are averaged with the segment widths and
areas to find the cell-face and cell-average values, respectively, for $I$ on
the underlying material spatial grid. These quantities are used to calculate
cell- and face-average angular moments of high-order transport solution
necessary to compute grid functions of the Eddington tensor in material
spatial cells. The cell-average Eddington tensor on the $i^{\text{th}}$
material grid cell is defined by

$\boldsymbol{\mathfrak{f}}_{i}=\frac{\sum_{m=1}^{M}\sum_{k,s\in
C(i,m)}\omega_{m}A_{k,s}(\boldsymbol{\Omega}_{m}\otimes\boldsymbol{\Omega}_{m})I_{k,s}}{\sum_{m=1}^{M}\sum_{k,s\in
C(i,m)}\omega_{m}A_{k,s}I_{k,s}}\,,$ (15)

where $\omega_{m}$ are the angular quadrature weights and $C(i,m)$ is the set
of characteristic segments for the direction $\boldsymbol{\Omega}_{m}$ which
trace over the $i^{\text{th}}$ material grid cell. The face-average Eddington
tensor on the $f^{\text{th}}$ cell face in the material grid is

$\boldsymbol{\mathfrak{f}}_{f}=\frac{\sum_{m=1}^{M}\sum_{k,s\in
F(f,m)}\omega_{m}w_{k}(\boldsymbol{\Omega}_{m}\otimes\boldsymbol{\Omega}_{m})I_{k,s}}{\sum_{m=1}^{M}\sum_{k,s\in
F(f,m)}\omega_{m}w_{k}I_{k,s}}\,,$ (16)

where $F(f,m)$ is the set of characteristic segments for the direction
$\boldsymbol{\Omega}_{m}$ whose upwind face intersects with the
$f^{\text{th}}$ material grid cell face.

The characteristic grid is constructed as follows. For each direction, an
initial set of characteristics are calculated such that no characteristic
traces over a vertex in the material spatial grid. The result is a ‘non-
uniform’ grid of characteristics for each direction with varying widths $w$.
We then enforce a maximum width $h_{\text{moc}}$. Any characteristics
calculated from the initial ray tracing procedure that have a width
$w>h_{\text{moc}}$ are split into 2 characteristics of width ${w}/{2}$.
Characteristics are continuously split into even halves by width until none
are wider than $h_{\text{moc}}$. Finally, each characteristic is formed into
segments length-wise whose boundaries are at the intersection point of a given
characteristic with a cell boundary in the material grid.

With this formulation, the TRT problem will be discretized with two unique
computational grids in space. This yields more flexibility in that both the
underlying material grid and the high-order characteristic grid can be refined
independently of one another. If the characteristic grid is continuously
refined on a static material grid, the closures provided by the BTE solution
on the material grid are expected to converge to some value. However the
accuracy of these closures (and the low-order solution) will still be limited
by the material grid. When allocating computational resources in a simulation,
it may be unclear which grid will benefit the most from extra refinement. We
perform convergence studies of each grid independently in Section 4 to
quantify these effects.

##  4 Numerical Results

Analysis is performed using the classical Fleck-Cummings (F-C) test problem in
2D Cartesian geometry [28]. The test domain is a homogeneous $6\times 6$ cm
slab with material defined by the spectral opacity
$\varkappa_{\nu}=\frac{27}{\nu^{3}}(1-e^{-\nu/T})$ and energy density
$\varepsilon=c_{v}T$ where $c_{v}=0.5917a_{R}(T^{\text{in}})^{3}$.
$T^{\text{in}}=1$ KeV is the temperature of radiation incoming to the domain
on the left boundary, with the other boundaries being at vacuum. The slab is
initially at a temperature of $T^{0}=1$ eV. The F-C test is simulated for
times $t\in[0,3\text{ns}]$. 150 uniform time steps are used $\Delta t=2\times
10^{-2}$ ns. 17 frequency groups are used and 144 discrete directions. The
Abu-Shumays angular quadrature set q461214 with 36 discrete directions per
quadrant is used [29]. All low-order equations are discretized with the
backward-Euler time integrator and a second-order finite volumes scheme. We
define a uniform, orthogonal material spatial grid whose cells are width
$h_{\text{mat}}$.

Figure 1: Iterations required with the MLQD algorithm to reach a convergence
criteria of $\epsilon=10^{-14}$ at each time step while
$h_{\text{moc}}=10^{-3}$ cm and $h_{\text{mat}}=1.2,0.6,0.3,0.015$ cm.

Figure 1 shows the number of high-order transport iterations required with the
MLQD algorithm at each time step to reach a very tight convergence criteria of
$\epsilon=10^{-14}$. Such a tight criteria is enforced to ensure no effects
are missed in the iterative behavior, and comparisons of the solution on
different spatial grid sizes can be made to high precision. Results are shown
for $h_{\text{moc}}=10^{-3}$ cm with several values for $h_{\text{mat}}$.
Iterations increase slightly with refinements in $h_{\text{mat}}$. For all
considered mesh widths the iterations converge rapidly. Here $h_{\text{moc}}$
is refined to a high level for all considered material grids. For
$h_{\text{mat}}=0.15$ cm, there will be more than $150$ characteristics traced
over every material grid cell for each discrete direction. In this way the
characteristic grid can be considered static for all considered
$h_{\text{mat}}$. We note that iteration counts per time step change
insignificantly with changes in $h_{\text{moc}}$.

Tables 3, 3 and 3 display the results of mesh refinement studies conducted for
both $h_{\text{mat}}$ and $h_{\text{moc}}$. Each refinement study is conducted
with one of the mesh widths held static. In Table 3, $h_{\text{moc}}$ is held
constant at $10^{-3}$ cm while $h_{\text{mat}}$ is refined in several steps.
In Table 3, $h_{\text{mat}}$ is held constant at $0.15$ cm while
$h_{\text{moc}}$ is refined in several steps. Table 3 holds $h_{\text{mat}}$
at 0.6 cm while $h_{\text{moc}}$ is refined. Here $\|\Delta
y_{h}\|=\|y_{h}-y_{2h}\|_{L_{2}}$ and $\rho_{h}^{y}$ = $\frac{\|\Delta
y_{2h}\|}{\|\Delta y_{h}\|}$. $\rho_{h}^{y}$ signifies the estimated
convergence rate of the variable $y$ with mesh refinement. For Table 3,
$h_{\text{moc}}=10^{-3}$ cm was chosen so that refinements in $h_{\text{mat}}$
would not cause implicit refinements in the characteristic grid. This allows
for analysis on the effect of strictly refining the material grid with an
effectively constant characteristic grid. $h_{\text{mat}}=0.15$ cm was chosen
as the material grid width in Table 3 to give a reasonably fine mesh on which
to study characteristic grid refinements. Similarly, $h_{\text{mat}}=0.6$ was
chosen as a suitably coarse material grid for Table 3.

Table 1: Convergence behavior for $T$ and $E$ with refinements in
$h_{\text{mat}}$ while $h_{\text{moc}}$ held constant at $10^{-3}$ cm

$h_{\text{mat}}$ | $\|\Delta T_{h_{\text{mat}}}\|$ | $\rho_{h_{\text{mat}}}^{T}$ | $\|\Delta E_{h_{\text{mat}}}\|$ | $\rho_{h_{\text{mat}}}^{E}$
---|---|---|---|---
0.6 | 19.1 | - | 2.13$\times 10^{-1}$ | -
0.3 | 9.70 | 1.96 | 1.10$\times 10^{-1}$ | 1.94
0.15 | 4.88 | 1.99 | 5.69$\times 10^{-2}$ | 1.97

Table 2: Convergence behavior for $T$ and $E$ with refinements in
$h_{\text{moc}}$ while $h_{\text{mat}}$ held constant at $0.15$ cm

$\frac{h_{\text{mat}}}{h_{\text{moc}}}$ | $h_{\text{moc}}$ | $\|\Delta T_{h_{\text{moc}}}\|$ | $\rho_{h_{\text{moc}}}^{T}$ | $\|\Delta E_{h_{\text{moc}}}\|$ | $\rho_{h_{\text{moc}}}^{E}$
---|---|---|---|---|---
8 | $2\times 10^{-2}$ | $8.01\times 10^{-5}$ | - | $8.71\times 10^{-7}$ | -
16 | $1\times 10^{-2}$ | $3.47\times 10^{-5}$ | 2.31 | $3.08\times 10^{-7}$ | 2.82
32 | $5\times 10^{-3}$ | $2.58\times 10^{-5}$ | 1.35 | $2.01\times 10^{-7}$ | 1.53
64 | $2.5\times 10^{-3}$ | $7.57\times 10^{-6}$ | 3.41 | $7.51\times 10^{-8}$ | 2.68

Table 3: Convergence behavior for $T$ and $E$ with refinements in
$h_{\text{moc}}$ while $h_{\text{mat}}$ held constant at $0.6$ cm

$\frac{h_{\text{mat}}}{h_{\text{moc}}}$ | $h_{\text{moc}}$ | $\|\Delta T_{h_{\text{moc}}}\|$ | $\rho_{h_{\text{moc}}}^{T}$ | $\|\Delta E_{h_{\text{moc}}}\|$ | $\rho_{h_{\text{moc}}}^{E}$
---|---|---|---|---|---
8 | $7.5\times 10^{-2}$ | $5.07\times 10^{-3}$ | - | $5.45\times 10^{-5}$ | -
16 | $3.75\times 10^{-2}$ | $1.96\times 10^{-3}$ | 2.59 | $2.59\times 10^{-5}$ | 2.10
32 | $1.875\times 10^{-2}$ | $1.04\times 10^{-3}$ | 1.88 | $1.05\times 10^{-5}$ | 2.48
64 | $9.375\times 10^{-3}$ | $2.93\times 10^{-4}$ | 3.55 | $3.04\times 10^{-6}$ | 3.44
128 | $4.6875\times 10^{-3}$ | $1.56\times 10^{-4}$ | 1.88 | $1.04\times 10^{-6}$ | 2.93
256 | $2.34375\times 10^{-3}$ | $4.35\times 10^{-5}$ | 3.59 | $2.99\times 10^{-7}$ | 3.48

The convergence rates for both $T$ and $E$ displayed for refinements in
$h_{\text{mat}}$ in Table 3 are very close to first order. The convergence
rates for refinements in $h_{\text{moc}}$ as shown in Tables 3 and 3 are not
as clear cut, taking on values between 1.3 and 3.6. Figures 2 and 3 display
plots of $\|\Delta T_{h_{\text{moc}}}\|$ and $\|\Delta E_{h_{\text{moc}}}\|$
in log-log format. The values shown in both Tables 3 and 3 are plotted
together, and separately alongside sample lines representative of first and
second order convergence rates. Visually, the solution converges at a rate
closer to first order when $h_{\text{mat}}=0.15$ cm, and closer to second
order while $h_{\text{mat}}=0.6$ cm.

Figure 2: Plots of $\|\Delta T_{h_{\text{moc}}}\|$ vs $h_{\text{moc}}$ from Tables 3 and 3 Top: $h_{\text{mat}}=0.6,\ 0.15$ cm; Bottom left: $h_{\text{mat}}=0.15$ with sample convergence rate lines; Bottom right: $h_{\text{mat}}=0.6$ with sample convergence rate lines.  Table 4: Relative differences in $T$ and $E$ from refinements in $h_{\text{moc}}$ while $h_{\text{mat}}$ held constant at $0.6$ cm $\frac{h_{\text{mat}}}{h_{\text{moc}}}$ | $h_{\text{moc}}$ | $\|\Delta T_{h_{\text{moc}}}\|_{\text{rel}}$ | $\|\Delta E_{h_{\text{moc}}}\|_{\text{rel}}$
---|---|---|---
8 | $7.5\times 10^{-2}$ | $9.75\times 10^{-5}$ | $1.53\times 10^{-4}$
16 | $3.75\times 10^{-2}$ | $3.76\times 10^{-5}$ | $7.31\times 10^{-5}$
32 | $1.875\times 10^{-2}$ | $2.01\times 10^{-5}$ | $2.95\times 10^{-5}$
64 | $9.375\times 10^{-3}$ | $5.64\times 10^{-6}$ | $8.57\times 10^{-6}$
128 | $4.6875\times 10^{-3}$ | $3.00\times 10^{-6}$ | $2.93\times 10^{-6}$
256 | $2.34375\times 10^{-3}$ | $8.36\times 10^{-7}$ | $8.43\times 10^{-7}$

The observed effects in $h_{\text{moc}}$ convergence behaviors could partially
be attributed to the fact that in constructing the characteristic mesh, only a
maximum width of characteristics is enforced. Depending on how the initial
grid is constructed, halving this maximum width criteria may not result in a
strict doubling of the number of characteristics. Furthermore, this effect can
change with the direction of motion $\boldsymbol{\Omega}$. At the moment,
there does not exist any rigorous theory which gives an expected rate of
convergence for the solution with refinements in the characteristic mesh. The
results given here show that the scheme converges but we do not see asymptotic
behavior for the rate of convergence. To fully understand the MOLC/RTS schemes
for the class of problems at hand, more detailed analysis will be required.

Figure 3: Plots of $\|\Delta E_{h_{\text{moc}}}\|$ vs $h_{\text{moc}}$ from Tables 3 and 3 Top: $h_{\text{mat}}=0.6,\ 0.15$ cm; Bottom left: $h_{\text{mat}}=0.15$ with sample convergence rate lines; Bottom right: $h_{\text{mat}}=0.6$ with sample convergence rate lines.  Table 5: Relative differences in $T$ and $E$ from refinements in $h_{\text{moc}}$ while $h_{\text{mat}}$ held constant at $0.15$ cm $\frac{h_{\text{mat}}}{h_{\text{moc}}}$ | $h_{\text{moc}}$ | $\|\Delta T_{h_{\text{moc}}}\|_{\text{rel}}$ | $\|\Delta E_{h_{\text{moc}}}\|_{\text{rel}}$
---|---|---|---
8 | $2\times 10^{-2}$ | $6.18\times 10^{-6}$ | $9.71\times 10^{-6}$
16 | $1\times 10^{-2}$ | $2.68\times 10^{-6}$ | $3.43\times 10^{-6}$
32 | $5\times 10^{-3}$ | $1.99\times 10^{-6}$ | $2.24\times 10^{-6}$
64 | $2.5\times 10^{-3}$ | $5.84\times 10^{-7}$ | $8.37\times 10^{-7}$

Tables 4 and 5 give the relative normed difference in solutions between
subsequent $h_{\text{moc}}$ grids. The relative differences are calculated as
$\|\Delta
y_{h}\|_{\text{rel}}={\|y_{h}-y_{2h}\|_{L_{2}}}\big{/}{\|y_{h}\|_{L_{2}}}$.
The relative error norms in both $T$ and $E$ are shown to decrease by roughly
an order of magnitude from $h_{\text{mat}}=0.6$ to $h_{\text{mat}}=0.15$ for
the same ratio of $\frac{h_{\text{mat}}}{h_{\text{moc}}}$. Instead if the two
tables are compared at similar values for $h_{\text{moc}}$ widths, the
relative difference values are similar to one another. In any case, the
magnitudes of both the absolute and relative error norms are small. Absolute
values are capped at $5\times 10^{-3}$ and relative values are capped at
$10^{-4}$. In comparison to the results in Table 3, the error norms are
significantly smaller for refinements in the characteristic grid.

When considering how to formulate the computational mesh for a simulation, it
is important to allocate resources such that the grid which limits solution
accuracy is refined to the highest level. The results presented here
demonstrate that, at least for the refinement levels examined, it is the
underlying material grid which limits solution accuracy. Using
$\frac{h_{\text{mat}}}{h_{\text{moc}}}=8$ for both $h_{\text{mat}}=0.6,\ 0.15$
cm is enough to ensure the characteristic grid is not the limiting grid in
accuracy. It is reasonable to extrapolate that
$\frac{h_{\text{mat}}}{h_{\text{moc}}}=8$ will continue to be satisfactory for
finer values of $h_{\text{mat}}$ given that the error norms in Tables 3 and 3
decrease with smaller $h_{\text{mat}}$ for the same
$\frac{h_{\text{mat}}}{h_{\text{moc}}}$.

##  5 Conclusions

In this paper the MLQD method with MOLC/RTS for the BTE is formulated for TRT
problems and analyzed. The method was tested on the 2D F-C test problem. The
analyzed method presents an independent discretization scheme for HEDP
simulations which allows for straightforward refinement or coarsening of the
discrete BTE solution separately from involved multiphysics equations. As the
BTE solution is refined, closures for the LOQD equations will become more
resolved on the given material grid. Iterations are shown to converge rapidly
at every time step for several refinements of the material grid. Independent
mesh refinement studies were performed for the underlying material spatial
grid and the characteristic grid. First order convergence of the solution for
$T$ and $E$ is observed for refinements in the material grid. Refinements in
the characteristic grid for the BTE do not give a clear convergence rate,
instead alternating between sub- and super-linear rates for $T$ and $E$. This
could be due to the fact that refinements in the characteristic grid do not
strictly double the number of characteristics. A more detailed analysis of
these effects should be partaken in the future to understand the convergence
behavior of the MOLC/RTS in these problems. The absolute difference between
solutions for subsequent mesh refinements on the material grid was shown to be
significantly larger than for refinements of the characteristic grid.

## Acknowledgements

Los Alamos Report LA-UR-22-32224. This research project was funded by the
Sandia National Laboratory, Light Speed Grand Challenge, LDRD, Strong Shock
Thrust. The work the first author (JMC) was also supported by the U.S.
Department of Energy through the Los Alamos National Laboratory. Los Alamos
National Laboratory is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of U.S. Department of Energy
(Contract No. 89233218CNA000001). The content of the information does not
necessarily reflect the position or the policy of the federal government, and
no official endorsement should be inferred.

## References

  * [1] Y. B. Zeldovich, Y. P. Razier, Physics of Shock Waves and High Temperature Hydrodynamic Phenomena, Academic, New York, 1966.
  * [2] R. P. Drake, High Energy Density Physics: Fundamentals, Inertial Fusion and Experimental Astrophysics, Springer, 2006.
  * [3] V. Vladimirov, Numerical solution of the kinetic equation for a sphere, Computational Mathematics 3 (1958) 3, (in Russian).
  * [4] V. Ya. Gol’din, Characteristic difference scheme for non-stationary kinetic equation, Soviet Mathematics Doklady 1 (1960) 902–906.
  * [5] K. Takeuchi, A numerical method for solving the neutron transport equation in finite cylindrical geometry, Journal of Nuclear Science and Technology 6 (1969) 446–473.
  * [6] J. R. Askew, A characteristic formulation of the neutron transport equation in complicated geometries, Tech. Rep. M1108, Atomic Energy Establishment (1972).
  * [7] H. D. Brough, C. T. Chudley, Characteristic ray solutions of the transport equation, in: J. Lewis, M. Becker (Eds.), Advances in Nuclear Science and Technology, Vol. 12, Plenum Press, 1980, pp. 1–30.
  * [8] R. E.Alcouffe, E.W. Larsen, A review of characteristic methods used to solve the linear transport equation, in: Proc. Int. Topl. Mtg. Advances in Mathematical Methods for the Solution of Nuclear Engineering Problems, Vol. 1, Munich, Germany, 1981, p. 3.
  * [9] J. R. Askew, M. J. Roth, WIMS-E: A scheme for neutronics calculations, Tech. Rep. AEEW-R-1315, United Kingdom Atomic Energy Authority (1982).
  * [10] M. Edenius, K. Ekberg, B. H. Forssén, D. Knott, CASM0-4, A Fuel Assembly Burnup Program, User’s manual, Studsvik/SOA-93/1, Studsvik of America, 1993.
  * [11] K. S. Smith, J. D. Rhodes, III, Full-core, 2-D, LWR core calculations with CASMO-4E, in: Proceedings of PHYSOR 2002, Seoul, Korea, 2002, p. 12 pp.
  * [12] R. Sanchez, A. Chetaine, A synthetic acceleration for a two dimensional characteristic method in unstructured meshes, Nuclear Science and Engineering 136 (2000) 122–139.
  * [13] M. Zika, M. Adams, Acceleration for long-characteristics assembly-level transport problems, Nuclear Science and Engineering 134 (2000) 135–158.
  * [14] T. Abel, B. Wandelt, Adaptive ray tracing for radiative transfer around point sources, Monthly Notices of the Royal Astronomical Society 330 (2002) L53–L56.
  * [15] T. Pandya, M. Adams, Method of long characteristics applied in space and time, in: Proc. of International Conference on Mathematics, Computational Methods & Reactor Physics (M&C 2009), Saratoga Springs, NY, 2009.
  * [16] L. Buntemeyer, R. Banerjeer, T. Peters r, M. Klassendr, R. Pudritz, Radiation hydrodynamics using characteristics on adaptive decomposed domains for massively parallel star formation simulations, New Astronomy 43 (2016) 49–69.
  * [17] V. Ya. Gol’din, A quasi-diffusion method of solving the kinetic equation, USSR Comp. Math. and Math. Phys. 4 (1964) 136–149.
  * [18] L. H. Auer, D. Mihalas, On the use of variable Eddington factors in non-LTE stellar atmospheres computations, Monthly Notices of the Royal Astronomical Society 149 (1970) 65–74.
  * [19] V. Ya. Gol’din, B. N. Chetverushkin, Methods of solving one-dimensional problems of radiation gas dynamics, USSR Comp. Math. and Math. Phys. 12 (1972) 177–189.
  * [20] K.-H. A. Winkler, M. L. Norman and D. Mihalas, Implicit adaptive-grid radiation hydrodynamics, in: Multiple Time Scales, Academic Press, 1985, pp. 145–184.
  * [21] V. Ya. Gol’din, D. A. Gol’dina, A. V. Kolpakov, A. V. Shilkov, Mathematical modeling of hydrodynamics processes with high-energy density radiation, Problems of Atomic Sci. & Eng.: Methods and Codes for Numerical Solution of Math. Physics Problems 2 (1986) 59–88, in Russian.
  * [22] D. Mihalas, B. Weibel-Mihalas, Foundation of Radiation Hydrodynamics, Oxford University Press, 1984.
  * [23] D. Y. Anistratov, E. N. Aristova, V. Y. Gol’din, A nonlinear method for solving problems of radiation transfer in a physical system, Mathematical Modeling 8 (1996) 3–28, in Russian.
  * [24] E. N. Aristova, V. Ya. Gol’din, A. V. Kolpakov, Multidimensional calculations of radiation transport by nonlinear quasi-diffusion method, in: Proc. of Int. Conf. on Math. and Comp., M&C 1999, Madrid, Spain, 1999, pp. 667–676.
  * [25] D. Y. Anistratov, Stability analysis of a multilevel quasidiffusion method for thermal radiative transfer problems, Journal of Computational Physics 376 (2019) 186–209.
  * [26] D. Y. Anistratov, V. Ya. Gol’din, Nonlinear methods for solving particle transport problems, Transport Theory and Statistical Physics 22 (1993) 42–77.
  * [27] P. Ghassemi, D. Y. Anistratov, Multilevel quasidiffusion method with mixed-order time discretization for multigroup thermal radiative transfer problems, Journal of Computational Physics 409 (2020) 109315.
  * [28] J. A. Fleck, J. D. Cummings, An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport, J. of Comp. Phys. 8 (1971) 313–342.
  * [29] L. K. Abu-Shumays, Angular quadratures for improved transport computations, Transport Theory & Statistical Physics 30 (2001) 169–204.


11institutetext: IDLab, Dept. Electronics & Information Systems, Ghent
University – imec, Belgium
11email: {x.sittminoo, gerald.haesendonck<EMAIL_ADDRESS>22institutetext: KULeuven, Dept. Computer Science – Leuven.AI – Flanders Make,
Belgium
22email<EMAIL_ADDRESS>
# RMLStreamer-SISO: an RDF stream generator from streaming heterogeneous data

Sitt Min Oo 11 0000-0001-9157-7507 Gerald Haesendonck 11 0000-0003-1605-3855
Ben De Meester 11 0000-0003-0248-0987 Anastasia Dimou 22 0000-0003-2138-7972

###### Abstract

Stream-reasoning query languages such as CQELS and C-SPARQL enable query
answering over RDF streams. Unfortunately, there currently is a lack of
efficient RDF stream generators to feed RDF stream reasoners. State-of-the-art
RDF stream generators are limited with regard to the velocity and volume of
streaming data they can handle. To efficiently generate RDF streams in a
scalable way, we extended the RMLStreamer to also generate RDF streams from
dynamic heterogeneous data streams. This paper introduces a scalable solution
that relies on a dynamic window approach to generate RDF streams with low
latency and high throughput from multiple heterogeneous data streams. Our
evaluation shows that our solution outperforms the state-of-the-art by
achieving millisecond latency (compared to seconds that state-of-the-art
solutions need), constant memory usage for all workloads, and sustainable
throughput of around 70,000 records/s (compared to 10,000 records/s that
state-of-the-art solutions take). This opens up the access to numerous data
streams for integration with the semantic web.

Resource type: Software

License: MIT License

URL: https://github.com/RMLio/RMLStreamer/releases/tag/v2.3.0

###### Keywords:

RML Stream processing Window Joins Knowledge graph generation

##  1 Introduction

An increasing portion of data are continuous in nature, e.g., sensor events,
user activities on a website, or financial trade events. This type of data is
known as data streams; sequences of unbounded tuples generated continuously in
different rates and volumes [3]. Due to the temporal nature of data streams,
low latency computation of analytical results is needed to timely react in
different use cases, e.g., fraud detection [9]. Thus, stream processing
engines must efficiently handle low latency computation of varying velocity
and volume.

On the one hand, different frameworks were proposed to handle data streams,
e.g., Flink, Spark or Storm [6, 26, 19]. On the other hand, RDF stream
processing (RSP) engines, e.g., CQELS and C-SPARQL [16, 1, 5], were widely
studied and perform high-throughput analysis of RDF streams with low memory
footprints [16]. Yet, these stream processing frameworks are not substantially
used in the domain of RDF graph generation from streaming data sources,
despite the demand of these mature RSP engines for more RDF streams.

Between data processing frameworks and stream processing engines, there are
tools to generate RDF streams from heterogeneous data streams (e.g. SPARQL-
Generate [17], RDFGen [21], TripleWave [18], Cefriel’s Chimera [22]). However,
some of these tools are inefficient when the data stream starts to scale in
terms of volume and velocity, such as TripleWave, and SPARQL-Generate. While
other tools are not open sourced nor suitable for the mapping of streaming
data, such as RDFGen, and Cefriel’s Chimera respectively. Overall, there are
no RDF stream generators that keep up with the needs of stream reasoning
engines while taking advantage of data processing frameworks to efficiently
produce RDF streams.

In this paper, we present the RMLStreamer-SISO, a parallel, vertically and
horizontally scalable stream processing engine to generate RDF streams from
heterogeneous data streams of any format (e.g. JSON, CSV, XML, etc.). We
extended previous preliminary work [13] of heterogeneous data stream mapping
solution: an open source implementation on top of Apache Flink [6], available
under MIT license, which generates high volume RDF data from high volume
heterogeneous data. RMLStreamer-SISO extends RMLStreamer to also support any
input data streams and export RDF streams (Stream-In-Stream-Out (SISO)).
RMLStreamer-SISO now supports a much larger part of the RML
specification111Implementation report of RML: https://rml.io/implementation-
report/, including all features of RML but relational databases.

The RMLStreamer-SISO outperforms the the state-of-the-art tools when handling
high velocity data stream, increasing the throughput it could handle while
maintaining low latency. The RMLStreamer-SISO achieves millisecond latency, as
opposed to seconds that state-of-the-art solutions need, constant memory usage
for all workloads, and sustainable throughput of around 70,000 records/s,
compared to 10,000 records/s that state-of-the-art solutions take.

Through the utilization of a low-latency tool like RMLStreamer-SISO, legacy
streaming systems could exploit the unique characteristics of real-life
streaming data, while enabling analysts to exploit the semantic reasoning
using knowledge graphs in real-time and have access to more reliable data.

The contributions presented in this paper are: (i) an algorithm to generate
the RDF streams from heterogeneous streaming data; (ii) its implementation,
the RMLStreamer-SISO, as an extension of RMLStreamer; and (iii) an evaluation
demonstrating that the RMLStreamer-SISO outperform the state-of-the-art. The
paper is structured as follows: Section 2 discusses related work, Section 3
the approach and its implementation, Section 4 the evaluation of RMLStreamer-
SISO against state-of-the-art, Section 5 the results of our evaluations, and
Section 7 concludes our work with possible future works.

##  2 Related Works

Streaming RDF mapping engines transform heterogeneous data streams to RDF data
streams. Several solutions exist in the literature for generating RDF from
persistent data sources [23, 14, 13, 2], but only few generate RDF from data
streams [21, 17, 18]. Although the implementations details are elaborated in
these works, their evaluations are designed without considering the different
data stream behaviours nor the resource contention between different
evaluation components.

TripleWave [18] generates RDF streams from streaming or static data sources
using R2RML mappings, and publishes them as RDF stream. However, the R2RML
mappings of TripleWave are invalid according to the specifications of R2RML
and it does not support joins. Although it is purported to support several
input sources, the user has to write the code to process the input data and
iterate over them before using the tool. This can result in poor performance
from improper implementation. Last, it is not designed to support distributed
parallel processing, resulting in limited scaling with data volume and
velocity.

RDF-Gen [21] generates static or streaming RDF data from static or streaming
data sources. A Data connector communicates with the data source, iterates
over its data entries, and converts every entry to a record of values. These
records are converted to RDF using a graph template: a listing of RDF-like
statements with variables bound to the record values coming from data
connectors. RDF-Gen generates RDF on a per record basis, theoretically
allowing a distributed parallel processing set-up. However, the current
implementation and documentation show no indication of a clustered setup nor
how to run it.

SPARQL-Generate [17] extends SPARQL 1.1 syntax to support mapping of
heterogeneous data to RDF data. SPARQL-Generate could be implemented on top of
any SPARQL query engine, and knowledge engineers with SPARQL experience could
use it with ease. The reference implementation of SPARQL-Generate222SPARQL-
Generate: https://github.com/sparql-generate/sparql-generate generates RDF
streams from data streams, even though it is not reported in the original
paper. Although joining data from multiple sources is supported, SPARQL-
Generate waits for one of the data streams to end first before consuming other
data sources to join the data. Thus, joins with unbounded streaming data
sources are not supported. The implementation is based on single machine setup
without scaling with data volume and velocity.

Cefriel’s Chimera [22] is an integration framework based on Apache Camel
333Apache Camel: https://camel.apache.org/ split into four “blocks” of
components to map heterogeneous data to RDF data: lifting block, data
enricher, inference enricher, and lowering block. Chimera aims to be modular
and allows each block to be replaced with custom implementations. The current
implementation uses a modified version of RMLMapper444RMLMapper:
https://github.com/RMLio/rmlmapper-java in the lifting block for data stream
processing. However, the whole RML mapping process is recreated with each
incoming message which could lead to high performance overhead in a highly
dynamic data stream environment.

##  3 Stream In - Stream Out (SISO)

We extend RMLStreamer’s architecture for generating RDF from persistent big
data sources [13] to also generate RDF streams from heterogeneous data streams
with high data velocity and volume, while keeping the latency low. The RDF
mapping language (RML) [10], a superset of R2RML, expresses customized mapping
from heterogeneous data sources to RDF datasets. We illustrate the concepts of
RML with the example RML document in Listing 2.

We break the process of generating RDF from a data stream into tasks and
subtasks (Figure 1). Each task or subtask is a stream processing operator
acting on an incoming data stream. They could be chained one after the other
to form a pipeline of operators and result in one or more outgoing data
streams. This approach introduces parallelism on both data and processing
level, enabling each data stream and operator to be processed and executed
respectively in parallel.

To illustrate RMLStreamer-SISO’s pipeline, we use the examples in Listing 1
and 2. The mapping document in Listing 2 is used to join and map JSON data
(Listing 1) from websocket streams to RDF with dynamic window join.

Figure 1: Workflow of RMLStreamer. Data flows from the Data Source at the top
through all the components pipeline to the Sink writer at the bottom. Listing
1: Data records from 2 data streams “Flow” & “Speed”.

⬇

1// data records from Speed stream

2{"speed":123.0,"time":"14:42:00","id":"lane1"}

3// data records from Flow stream

4{"flow":1680,"time":"14:42:00","id":"lane1"}

Listing 2: Example RML Mapping file to generate streaming RDF from the
streaming heterogeneous data of Listing 1.

⬇

1# prefix definitions omitted

2_:ws_source_ndwSpeed a td:Thing ;

3 td:hasPropertyAffordance [ td:hasForm [

4 hctl:hasTarget "ws://data-streamer:9001" ; # URL and content type

5 hctl:forContentType "application/json" ; # Data format

6 hctl:hasOperationType "readproperty" ] ] . # Read only

7_:ws_source_ndwFlow a td:Thing;

8 td:hasPropertyAffordance [ td:hasForm [

9 hctl:hasTarget "ws://data-streamer:9000" ;

10 hctl:forContentType "application/json" ;

11 hctl:hasOperationType "readproperty" ] ] .

12<JoinConfigMap> a rmls:JoinConfigMap ;

13 rmls:joinType rmls:TumblingJoin . # Trigger/eviction type

14<NDWSpeedMap> a rr:TriplesMap ;

15 rml:logicalSource [ # Describes data source

16 rml:source _:ws_source_ndwSpeed ;

17 rml:referenceFormulation ql:JSONPath ; # JSONPath iterator

18 rml:iterator "$" ] ; # Iterates the data as JSON root object

19 rr:subjectMap [ # Generation of the subject IRI

20 rr:template "speed={speed}&time={time}" ] ;

21 rr:predicateObjectMap [ # Describes how predicate and object are generated

22 rr:predicate <http://example.com/laneFlow> ;

23 rr:objectMap [

24 rr:parentTriplesMap <NDWFlowMap> ; # TripleMap to be joined with

25 rmls:joinConfig <JoinConfigMap> ; # Configuration of join window

26 rmls:windowType rmls:TumblingWindow ; # Type of join window

27 rr:joinCondition [ # Attributes on which the data records are joined

28 rr:child "id" ; rr:parent "id" ; ] ] ] .

29<NDWFlowMap> a rr:TriplesMap ;

30 rml:logicalSource [

31 rml:source _:ws_source_ndwFlow ;

32 rml:referenceFormulation ql:JSONPath ;

33 rml:iterator "$" ] ;

34 rr:subjectMap [ rr:template "flow={flow}&time={time}" ] .

###  3.1 RDF stream generation workflow

Our approach consists of a workflow with four tasks (see Figure 1):

#### Ingestion

The ingestion task captures data streams and prepares the data records for the
mapping task. Each data stream triggers one ingestion task that can run in
parallel with the other ingestion tasks spawned by the other data streams. The
ingestion task can be divided in three subtasks:

  1. 1.

Data source connector (Figure 1, (b)): This subtask is responsible for
connecting to a (streaming) data source (a). It reads data records from the
source and passes these records (c) on to the stream partitioner.

  2. 2.

Stream partitioner (d): The stream of data records is optionally partitioned
in disjoint partitions to be fed to the next subtask. The partitioning depends
on the order’s maintenance. If the exact order of the incoming data records is
not important to be maintained, then these records can be distributed evenly
among multiple instances of the next subtask, increasing parallelism. If the
order of generating RDF statements needs to correspond with the order of the
incoming data records, then the stream is not distributed at this stage.

  3. 3.

Item generator (e): One data record can lead to zero or more RDF statements.
This subtask splits a data record in zero or more items of internal
representation called data items (f), according to the logical iterators
defined in the mapping document, before the actual mapping task takes place.
Using the sample data and the mapping document from Listing 1 and 2
respectively, this subtask will use the logical iterator ‘$’, a
JSONPath555JSONPath documentation:
https://goessner.net/articles/JsonPath/index.html, to generate data items from
each data record shown in Listing 1. In this case, the logical iterator is the
JSON root object, so the data item is the same as the incoming data records.
Otherwise, if the data record contains a list of sub-records, and the logical
iterator is specified over the list (e.g., $.list[*]), each of these sub-
records are turned into data items.

#### Pre-mapping (optional)

Before the data items are mapped to RDF, the data items may be processed with
custom data transformations defined with FnO [8], or the window operators,
such as joins, aggregates, and reduce. The FnO functions could be as simple as
changing letters to uppercase or as complex as the window joins. This stage is
optional and omitted if the RML document does not define pre-mapping
functions. The pre-mapping task (g) is right before the mapping task since the
data fields requiring preprocessing can be more than the data fields needed
for mapping to RDF data. For example, with the given inputs and mapping
document (Listing 2), the data items (Listing 1) from the two input streams,
“Flow” and “Speed”, are first buffered inside a window, and then joined based
on their internalId value. Data records, having the same value for the “id”,
are joined pairwise. If windows joins were implemented after the mapping
stage, the verbosity of RDF would substantially increase the network
bandwidth. More, to fully map the data before joining, RMLStreamer-SISO needs
to know all attributes present in the raw data records which would be
infeasible.

To support joining with windows, RML was extended. New vocabulary terms were
defined to support windowing operations with RML. We defined two new
properties: rmls:windowType to provide the type of window to be used when
joining and rmls:joinConfig when joining the Child and Parent Triple Map to
define how the trigger, and eviction are fired inside the window.

Section 3.2 details the dynamic windowing algorithm and Section 3.3 elaborates
on the design choice and windows’ implementation for RMLStreamer-SISO.

#### Mapping

RDF statements are generated from data items coming from the ingestion task
and the pre-mapping task.

  1. 1.

Statement generator (h): Each data item leads to one or more RDF statements in
this sub task. Each statement is generated in parallel as an abstract RDF
statement (i) which could be fed to the next subtask for serialization.

  2. 2.

RDF serializer (j): The abstract RDF statements are serialized into various
RDF serialisations based on the configuration given to the RMLStreamer.

#### Combination

This task brings back together all streams of RDF statements (l) into one
final RDF stream which will be written out using the sink writer (m).

###  3.2 Heterogeneous data streams join in RDF streams

Supporting joins in RMLStreamer-SISO and any streaming RDF generator, is not
trivial as windowing techniques are required for unbounded and unsynchronized
streaming data. Unlike batch processing where data is bounded, processing
whole data streams in memory is unsustainable due to the continuous and
infinite characteristics of streaming data. Therefore, stream processing
engines use buffers called windows to hold the most recent stream of records
in memory. The windows’ lifetime is measured in terms of time interval, thus,
the window interval determines the size of the window and their operation
behaviour is defined by the trigger, and the eviction events [12]. A trigger
event occurs when an operator is executed to process the data records inside
the window interval. An eviction event occurs when the window evicts the data
records inside its buffer.

Figure 2: Behaviour of the dynamic window under high, and low stream
velocities. The cogwheels are the trigger events representing the moment when
the data records are processed. In this figure, the trigger events are fired
with every new data record, and only when there is at least one data record
from each data stream.

We opted for an eager trigger implementation to lower the latency of
RMLStreamer’s responses for the windowed joins’ implementation. The joined
results are emitted as soon as possible without waiting for the eviction event
to occur. We designed a dynamic window which adapts its window intervals
according to the velocity of the incoming data streams. Adaptive windowing
[27] was studied in the context of batch stream processing with a positive
impact on the stream processing job’s performance: lower latency, and higher
throughput. We opted for a simple cost metric based on the data records’
number to keep the memory and latency low in a real-time stream processing
environment where the time constraint is more stringent.

The algorithm is inspired by the additive-increase, and multiplicative
decrease algorithm of TCP congestion control [7]. Figure 2 shows the high
level behaviour of our dynamic window for the two different stream velocities.
When the data stream velocity is high, the window size shrinks to process the
data records as fast as possible, keeping the latency low and throughput high.
When the data stream velocity is low, the size of the window grows to wait for
more data records and process them. This ensures that the window do not miss
the records due to short window size. We elaborate the details of the
algorithm below. For each window, the following configuration parameters are
provided:

  1. 1.

$|W|$: The window interval

  2. 2.

$\epsilon_{u}$ & $\epsilon_{l}$: Upper and lower threshold limit for total
cost metric

  3. 3.

$U$ & $L$: Upper and lower limit for the window interval

  4. 4.

$Limit(List_{P})$ & $Limit(List_{C})$: Upper limit size for parent and child
stream

Data:
$|W|,\epsilon_{u},\epsilon_{l},U,L,Limit(List_{P}),Limit(List_{C}),S_{P},S_{C}$

1 $cost(List_{P})=|S_{P}|/Limit(List_{P})$

2$cost(List_{C})=|S_{C}|/Limit(List_{C})$

3total cost $m=cost(List_{P})+cost(List_{C})$

// adapts window size based on cost

4 if  _$m >\epsilon_{u}$_ then

5 $|W|=|W|/2$

6 $Limit(List_{P})=Limit(List_{P})*cost(List_{P})*1.5$

7 $Limit(List_{C})=Limit(List_{C})*cost(List_{C})*1.5$

8else if  _$m <\epsilon_{l}$_ then

9 $|W|=|W|*1.1$

10 $Limit(List_{P})=Limit(List_{P})*cost(List_{P})*1.5$

11 $Limit(List_{C})=Limit(List_{C})*cost(List_{C})*1.5$

12

13clean both $List_{C}$ and $List_{P}$

clip $|W|$ in the range of $[L,U]$

Algorithm 1 Dynamic window $onEviction$ routine

Since we implement the join operator with eager execution, the trigger event
is fired when the current record $r_{c}$ arrives inside the window. We denote
the current window as $W$ with interval size $|W|$. The streams are denoted as
$S_{p}$ and $S_{c}$ with the corresponding states $List_{p}$ and $List_{c}$,
for the parent and child stream respectively (the parent and child stream
follows the RML specification for joining triples maps). The states contain
the records from their respective streams inside the window with for example
$|S_{p}|$ denoting the number of records from $S_{p}$. $List_{p}$ and
$List_{c}$ are only used in cost calculation to determine if the window
interval needs to be changed; they do not limit the amount of records that
could be buffered inside the window.

At each eviction trigger, we calculate the cost for each list states
$List_{p}$ and $List_{c}$. For example, the cost for
$cost(List_{P})=|S_{p}|/Limit(List_{P})$. The total cost is
$m=cost(List_{p})+cost(List_{c})$ and it is checked against the thresholds
$\epsilon_{l}$ and $\epsilon_{u}$ to adjust the window interval accordingly.
We assume the stable zone to be achieved if the total cost fulfils the
predicate: $\epsilon_{l}\leq m\leq\epsilon_{u}$. Algorithm 1 shows the pseudo-
code for the eviction algorithm we just elaborated.

###  3.3 Implementation

RMLStreamer-SISO is released as version 2.3 of RMLStreamer to utilize Flink’s
parallelism for horizontal scaling (via distributed processing in a network
and vertical scaling (via multi-threaded execution of tasks). The update
brings the windowing support for joining multiple data streams, the dynamic
windowing algorithm, and FnO [8] as an extension point for joins execution.
The code and usage instructions for RMLStreamer-SISO are available online at
the Github repository: https://github.com/RMLio/RMLStreamer.

Windowing support is implemented through the use of Flink’s windowing
API666Window: https://nightlies.apache.org/flink/flink-docs-
release-1.14/docs/dev/datastream/operators/windows/ for common types of
window, e.g., Tumbling Window. We implemented the KeyedCoProcessFunction
provided by Flink’s low-level stream processing API to manage the different
states required for the algorithm (Algorithm 1) of the dynamic window. We
implemented the dynamic windowed join before the mapping stage, to group input
streams and reduce network bandwidth usage. The generated RDF stream could be
windowed by the RDF stream processing engines consuming the output.

Currently, FnO functions jar files have to be compiled together as part of the
RMLStreamer-SISO jar. Examples on github777RMLStreamer-SISO:
https://github.com/RMLio/RMLStreamer show the working of RMLStreamer-SISO with
TCP data stream. We also provide an extensive documentation on RMLStreamer-
SISO in a containerized environment with docker888Docker: https://docker.com.

##  4 Evaluation

An extensive evaluation was conducted focused on variable data stream
velocity, volume and variety of data formats to emulate the real-life
workloads as close as possible. The code for the evaluation is available on
github999Benchmark: https://github.com/s-minoo/rmlstreamer-benchmark-rust.
Since RMLStreamer-SISO is situated between traditional stream processing and
RSP, state-of-the-art approaches for benchmarking in these domains are
combined: architectural design of RSPLab [25], workload design of Open Stream
Processing Benchmark [11], and measurement strategies of Karimov et. al [15].

We compare the RMLStreamer-SISO with the state-of-the-art streaming RDF
generator, SPARQL-Generate, which is actively maintained, used and supports
the same features as RMLStreamer-SISO. The other tools were not considered for
different reasons: TripleWave requires a custom implementation to process each
data stream and feed it in TripleWave which means it cannot be used as-it-is.
More, TripleWave is meant purely for feeding RDF streams to RDF stream
processing engines without performing joins, therefore it would have been an
unfair comparison both in terms of features and scope. RDF-Gen’s source code
is not available, but only a jar is available without any instructions to run
it. Both TripleWave and RDF-Gen are also not actively maintained. Finally,
Cefriel’s Chimera restarts the RDF mapping engine with every data record,
which means that the processing of the input and mapping is not performed in a
true streaming manner; the comparison would not be meaningful.

#### Data source

The input data used in the evaluation comes from time annotated traffic sensor
data from the Netherlands provided by NDW (Nationale Databank
Wegverkeersgegevens)101010NDW: https://www.ndw.nu/, and also used by Van
Dongen et al [11]. It contains around 68,000 rows of CSV data with two
different measurements across different lanes on a highway: number of cars
(flow), and their average speed. The two measurements are streamed through a
websocket data streaming server.

#### Metrics

Stream processing frameworks are typically evaluated using two main metrics:
latency and throughput [15]. Latency can be further distinguished into two
types: processing-time latency, and event-time latency. Processing time
latency is the interval between the data record’s arrival time at the input
and the emission time at the output of the streaming engine [15]. Event-time
latency is the interval between the creation time, and the emission time at
the streaming engine’s output, of the data record [15]. Latency measurement
requires to consider the effect of coordinated occlusion, where the queueing
time, a part of the event-time, is ignored [15]. Therefore, we consider event-
time latency as our latency measurement to take the effect of coordinated
occlusion in consideration.

For our evaluation, we considered the event-time latency of each record, the
throughput as number of consumed records per second, the memory and CPU usage
of the engine’s docker container. The measurements are captured on a machine
separate from the host machine of the System Under Test (SUT), where memory
and CPU usage are measured using cAdvisor111111cAdvisor:
https://github.com/google/cadvisor. By treating the SUTs as a blackbox, we
ensure that the measurement of the metrics incurs no performance penalty nor
resource contention with the SUTs during the evaluation.

Figure 3: Benchmark architecture to evaluate the different engines, inspired
by RSPLab.

#### Evaluation set up

The architectural design is a modification of RSPLab with a custom data
streaming component (Figure 3). It consists of three components: a) the data
streamer, b) the system under test, and c) the monitoring system.

With the proposed architecture where each components is isolated, we aim to
reduce the influence of the benchmark components on the engine during the
evaluation process. The modularity of the setup also increases the flexibility
of configuring the evaluation environment with minimal changes for the
engines.

#### Workload Design

To evaluate the performance of the engines under different data
characteristics and processing scenarios, we devise three different workloads:
(i) throughput measurement, (ii) periodic burst, and (iii) scalability
measurement. As SPARQL-Generate is unable to join unbounded streaming data (it
expects data streams with an end, Section 2), we evaluated the two workloads
(throughput measurement and periodic burst) without joining functionality to
compare.

  * •

throughput measurement: the data stream throughput is constant and steadily
increases with each run to determine the engine’s sustainable throughput [15].
CPU, latency and memory usage are measured.

  * •

periodic burst: a burst of data records is emitted periodically to mimic
fluctuations in data streams; CPU, memory, latency and throughput are
measured.

  * •

scalability measurement: RMLStreamer-SISO is evaluated in two modes:
centralised mode without parallelism and distributed mode with parallelizable
data to measure the impact of parallelism on its scalability. In both modes,
data from two input streams are joined and latency is measured.

#### System specifications

We ran the evaluation on a single machine with multiple docker containers to
emulate the communication between the data streaming source and the mapping
engine in a streaming network environment. The machine has Intel i7 CPU with 8
cores at 4.8Ghz, 16GB RAM, and 200GB hard disk space. The data streamer and
the monitoring system docker containers (Figure 3) have access to 4 of the
cores, and the SUT docker container has access to the leftover 4 cores. This
prevents CPU resource contention between the SUT and the other components used
for running the evaluation.

To evaluate horizontal scaling, the data streamer component is replaced with
Apache Kafka to support parallel ingestion of data streams by RMLStreamer-
SISO. Apache Kafka is configured with default settings and the data (Section
4) is streamed into two topics121212Kafka topics:
https://developer.confluent.io/learn-kafka/apache-kafka/topics/; “ndwFlow”,
and “ndwSpeed” containing the records about the number, and the average speed
of the cars respectively.

##  5 Results

In this section, we discuss the results of our evaluation using different
workloads.

#### Workload for throughput measurement

For the throughput measurement workload, we ran the evaluation multiple times
with increasing input data throughput for each run to evaluate the sustainable
throughput of the SUTs.

In the first few runs of the evaluation, the RMLStreamer fared a bit worse
than SPARQL-Generate in all three measurements. This is due to the overhead of
having a distributed task manager for executing, and managing the different
tasks and subtasks of mapping heterogeneous data (Figure 1). However, when the
throughput starts increasing beyond 10,000 records per second, RMLStreamer-
SISO outperforms SPARQL-Generate in terms of latency and memory usage.

Figure 4: SUTs performance under different data stream velocity for
sustainable throughput measurement. The last run for SPARQL-Generate was
omitted because it took more than 1 hour instead of the expected 30 minutes to
process the whole data stream.

Compared to RMLStreamer-SISO, SPARQL-Generate became unsustainable when the
throughput of the input data streams passes 10,000 records per second with 20s
latency (Figure 4). To the contrary, RMLStreamer-SISO has a consistent low
latency of 1 second for all runs of the workload having 100x magnitude lower
latency than SPARQL-Generate in later runs.

Regarding CPU usage, RMLStreamer-SISO has on average 20% more CPU usage for
the overhead of Apache Flink managing the distributed tasks.

In terms of memory usage, RMLStreamer-SISO uses significantly lower memory
than SPARQL-Generate at around 900MB compared to 3GB by SPARQL-Generate. Based
on the previous observations, we conclude that RMLStreamer-SISO outperforms
SPARQL-Generate at higher throughput with lower latency and memory usage. Even
though, RMLStreamer-SISO’s CPU usage is around 30% higher than SPARQL-Generate
in the last run, it effectively copes with the increase in data stream
velocity to maintain low latency processing.

#### Workload for periodic burst

The periodic burst workload studies the adaptability of the engine to the
recurring sudden burst of data stream. We used the measurements from the last
minute of the evaluation, when the engines are stable without warm-up
overheads, to better visualize their performance during the periodic burst of
data (Figure 5). In Figure 5, we see a periodic increase, and drop in the
throughput metrics measurements, which is an expected behaviour in the engines
when consuming a data stream input with periodic burst of data. Every 10
seconds we see a burst of around 35,000 messages. Both engines behave as
expected for the throughput metrics measurement.

The spikes for latency measurement of SPARQL-Generate (Figure 5) have a wider
base than those of RMLStreamer-SISO. This indicates that SPARQL-Generate takes
a longer time to recover from processing periodic workload than RMLStreamer-
SISO by a few seconds. RMLStreamer-SISO’s peak latency is around 500ms whereas
SPARQL-Generate has a peak latency of around $3.5s$. Although RMLStreamer-SISO
uses more CPU than SPARQL-Generate to process data burst, it adapts to the
sudden burst of data and recover more quickly than SPARQL-Generate. The
latency of RMLStreamer-SISO is also 7 times lower than SPARQL-Generate due to
the record-based processing capabilities. We conclude that RMLStreamer-SISO is
better adapted to workloads with periodic burst of data with faster recovery
period, lower latency and memory usage while maintaining the same throughput
capabilities as SPARQL-Generate.

Figure 5: Performance of SUTs in the last one minute of the periodic burst
workload evaluation. A part of the throughput graph is blurred to give more
clarity to the relationship between the trends in latency and throughput of
the engines.

#### Workload for scalability measurement

We evaluated the RMLStreamer-SISO’s capability to join two data streams with a
constant throughput of around 17000 messages per second. CPU and memory usage
of both modes of RMLStreamer-SISO are similar throughout the evaluation.
However, despite the similar performance in terms of CPU and memory usage,
parallelized mode fared significantly better in terms of the latency metric
than unparallelized mode. Unparallelized mode has a median latency of around
50000ms whereas parallelized mode has a median latency of around 57ms. This is
around 1000x lower in terms of the median latency. Moreover, the minimum
latency of parallelized RMLStreamer-SISO at 8ms is 10,000x lower than the
minimum latency of unparallelized RMLStreamer-SISO at 13653ms. The latency is
kept low with high parallelization due to the effective distribution of the
workload amongst the different parallelized tasks by the underlying DSP engine
(Apache Flink). We conclude that RMLStreamer scales extremely well with
significantly better performance in terms of latency if configured to be
executed in a distributed mode.

##  6 Use Cases

RMLStreamer-SISO has seen uptake in multiple projects – covering different use
cases in different architectures – to process streaming data and generate RDF
streams. Largest validation was in research and development (R&D) projects
between imec and Flemish companies such as DyVerSIFy on streaming data
analysis and visualisation [24, 20], together with Televic Rail on IoT data,
DAIQUIRI131313https://www.imec-int.com/en/what-we-offer/research-
portfolio/daiquiri together with VRT on sport sensor data, and ESSENCE and
H2020 project MOS2S141414https://innovatie.vrt.be/project/essence,
https://itea4.org/project/mos2s.html on media data. Other projects include
DiSSeCt151515https://smit.vub.ac.be/project/dissect on health data and
transport data [4]. The variety of use cases shows that the resource is
suitable for solving the task at hand and also applicable to a multitude of
use cases for society in general. Applications – within the knowledge graph
construction problem domain – are varied, i.e., processing a large amount of
low-frequency sensor data, a small amount of high-frequency sensor data, and
large data sets combined with streaming data, processing Kafka streams, MQTT,
Socket.io, and TCP streams. Beyond Belgium, RMLStreamer has received attention
by the Institute of Data Science, proposed as part of RDF graph generation
tutorials such as those by STIInnsbruck in Austria, and services such as
Data2Services161616https://maastrichtu-ids.github.io/best-
practices/blog/2021/03/18/build-a-
kg/,https://stiinnsbruck.github.io/lkgt/,https://d2s.semanticscience.org/docs/d2s-rml/
by the Institute of Data Science in Maastricht in the Netherlands.

##  7 Conclusion and Future Work

In this paper, we present RMLStreamer-SISO, a highly scalable solution to
seamlessly generate RDF streams thanks to its dynamic window algorithm which
adapts its window size to handle the dynamic characteristics of the data
stream. This way, RMLStreamer-SISO enables low latency and high throughput
mapping of heterogeneous data to RDF data. We showed that our solution scales
better than the state-of-the-art in terms of latency, memory, and throughput.
It is the only RDF stream generator which joins unbounded data streams and
scale horizontally and vertically, enabling RDF streams generation from
heterogeneous data streams which was not possible so far. Given it is open
source and already widely used in different use cases involving not only
academia but also industry, as shown in our use cases, it is expected that the
community that grew around it will further grow and contribute at its
maintenance, while its extensive documentation and tutorials allow for easy
reuse171717Example of tutorial for use with docker technology,
https://github.com/RMLio/RMLStreamer/tree/development/docker. The RML
extensions will be further discussed with the W3C community group on knowledge
graph construction and eventually will be incorporated to the revised RML
specification.

RMLStreamer-SISO increases the availability of RDF streams following the high
availability of data streams. Using a low-latency tool like RMLStreamer-SISO,
legacy streaming systems could exploit the unique characteristics of real-life
streaming data, while enabling analysts to exploit the semantic reasoning
using knowledge graphs in real-time. This way, we enabled access to more data
which should impact the further improvements of RSP engines and other semantic
web technologies on top of RDF streams which were not possible so far.

#### Resource Availability Statement:

Source code for RMLStreamer-SISO is available at
https://github.com/RMLio/RMLStreamer. The source code for the benchmark is
available at https://github.com/s-minoo/rmlstreamer-benchmark-rust. The
dataset used for the benchmark is available at
https://github.com/Klarrio/open-stream-processing-benchmark/tree/master/data-
stream-generator.

## References

  * [1] Barbieri, D.F., Braga, D., Ceri, S., Della Valle, E., Grossniklaus, M.: C-sparql: Sparql for continuous querying. In: Proceedings of the 18th International Conference on World Wide Web. p. 1061–1062. WWW ’09, Association for Computing Machinery, New York, NY, USA (2009). https://doi.org/10.1145/1526709.1526856
  * [2] Belcao, M., Falzone, E., Bionda, E., Valle, E.D.: Chimera: A bridge between big data analytics and semantic technologies. In: Hotho, A., Blomqvist, E., Dietze, S., Fokoue, A., Ding, Y., Barnaghi, P., Haller, A., Dragoni, M., Alani, H. (eds.) The Semantic Web – ISWC 2021. pp. 463–479. Springer International Publishing, Cham (2021)
  * [3] Botan, I., Derakhshan, R., Dindar, N., Haas, L., Miller, R.J., Tatbul, N.: Secret: A model for analysis of the execution semantics of stream processing systems. Proc. VLDB Endow. 3(1–2), 232–243 (sep 2010). https://doi.org/10.14778/1920841.1920874
  * [4] Brouwer, M.D., Bonte, P., Arndt, D., Sande, M.V., Heyvaert, P., Dimou, A., Verborgh, R., Turck, F.D., Ongenae, F.: Distributed continuous home care provisioning through personalized monitoring & treatment planning. In: Companion Proceedings of the Web Conference 2020. ACM (Apr 2020). https://doi.org/10.1145/3366424.3383528
  * [5] Calbimonte, J.P., Corcho, O., Gray, A.J.G.: Enabling ontology-based access to streaming data sources. In: Patel-Schneider, P.F., Pan, Y., Hitzler, P., Mika, P., Zhang, L., Pan, J.Z., Horrocks, I., Glimm, B. (eds.) The Semantic Web – ISWC 2010. pp. 96–111. Springer Berlin Heidelberg, Berlin, Heidelberg (2010)
  * [6] Carbone, P., Katsifodimos, A., Ewen, S., Markl, V., Haridi, S., Tzoumas, K.: Apache flink™: Stream and batch processing in a single engine. IEEE Data Eng. Bull. 38, 28–38 (2015)
  * [7] Chiu, D.M., Jain, R.: Analysis of the increase and decrease algorithms for congestion avoidance in computer networks. Computer Networks and ISDN systems 17(1), 1–14 (1989)
  * [8] De Meester, B., Dimou, A., Verborgh, R., Mannens, E.: An ontology to semantically declare and describe functions. In: Sack, H., Rizzo, G., Steinmetz, N., Mladenić, D., Auer, S., Lange, C. (eds.) The Semantic Web. pp. 46–49. Springer International Publishing, Cham (2016)
  * [9] Dias de Assunção, M., da Silva Veith, A., Buyya, R.: Distributed data stream processing and edge computing: A survey on resource elasticity and future directions. Journal of Network and Computer Applications 103, 1–17 (2018). https://doi.org/10.1016/j.jnca.2017.12.001, https://www.sciencedirect.com/science/article/pii/S1084804517303971
  * [10] Dimou, A., Vander Sande, M., Colpaert, P., Verborgh, R., Mannens, E., Van de Walle, R.: Rml: A generic language for integrated rdf mappings of heterogeneous data. vol. 1184 (04 2014)
  * [11] van Dongen, G., Van den Poel, D.: Evaluation of stream processing frameworks. IEEE Transactions on Parallel and Distributed Systems 31(8), 1845–1858 (2020). https://doi.org/10.1109/TPDS.2020.2978480
  * [12] Gedik, B.: Generic windowing support for extensible stream processing systems. Softw. Pract. Exper. 44(9), 1105–1128 (sep 2014). https://doi.org/10.1002/spe.2194
  * [13] Haesendonck, G., Maroy, W., Heyvaert, P., Verborgh, R., Dimou, A.: Parallel rdf generation from heterogeneous big data. In: Proceedings of the International Workshop on Semantic Big Data. SBD ’19, Association for Computing Machinery, New York, NY, USA (2019). https://doi.org/10.1145/3323878.3325802
  * [14] Iglesias, E., Jozashoori, S., Chaves-Fraga, D., Collarana, D., Vidal, M.E.: Sdm-rdfizer. Proceedings of the 29th ACM International Conference on Information & Knowledge Management (Oct 2020). https://doi.org/10.1145/3340531.3412881
  * [15] Karimov, J., Rabl, T., Katsifodimos, A., Samarev, R., Heiskanen, H., Markl, V.: Benchmarking distributed stream data processing systems. 2018 IEEE 34th International Conference on Data Engineering (ICDE) (Apr 2018). https://doi.org/10.1109/icde.2018.00169
  * [16] Le Phuoc, D., Dao-Tran, M., Le Tuan, A., Duc, M.N., Hauswirth, M.: Rdf stream processing with cqels framework for real-time analysis. In: Proceedings of the 9th ACM International Conference on Distributed Event-Based Systems. p. 285–292. DEBS ’15, Association for Computing Machinery, New York, NY, USA (2015). https://doi.org/10.1145/2675743.2772586
  * [17] Lefrançois, M., Zimmermann, A., Bakerally, N.: A sparql extension for generating rdf from heterogeneous formats. In: Blomqvist, E., Maynard, D., Gangemi, A., Hoekstra, R., Hitzler, P., Hartig, O. (eds.) The Semantic Web. pp. 35–50. Springer International Publishing, Cham (2017)
  * [18] Mauri, A., Calbimonte, J.P., Dell’Aglio, D., Balduini, M., Brambilla, M., Della Valle, E., Aberer, K.: Triplewave: Spreading rdf streams on the web. In: Groth, P., Simperl, E., Gray, A., Sabou, M., Krötzsch, M., Lecue, F., Flöck, F., Gil, Y. (eds.) The Semantic Web – ISWC 2016. pp. 140–149. Springer International Publishing, Cham (2016)
  * [19] N.A: Apache storm, https://storm.apache.org/
  * [20] Paepe, D.D., Hautte, S.V., Steenwinckel, B., Moens, P., Vaneessen, J., Vandekerckhove, S., Volckaert, B., Ongenae, F., Hoecke, S.V.: A complete software stack for IoT time-series analysis that combines semantics and machine learning—lessons learned from the dyversify project. Applied Sciences 11(24), 11932 (dec 2021). https://doi.org/10.3390/app112411932
  * [21] Santipantakis, G.M., Kotis, K.I., Vouros, G.A., Doulkeridis, C.: Rdf-gen: Generating rdf from streaming and archival data. In: Proceedings of the 8th International Conference on Web Intelligence, Mining and Semantics. WIMS ’18, Association for Computing Machinery, New York, NY, USA (2018). https://doi.org/10.1145/3227609.3227658
  * [22] Scrocca, M., Comerio, M., Carenini, A., Celino, I.: Turning transport data to comply with eu standards while enabling a multimodal transport knowledge graph. The Semantic Web – ISWC 2020 p. 411–429 (2020). https://doi.org/10.1007/978-3-030-62466-8_26
  * [23] Simsek, U., Kärle, E., Fensel, D.A.: Rocketrml - a nodejs implementation of a use case specific rml mapper. ArXiv abs/1903.04969 (2019). https://doi.org/10.48550/ARXIV.1903.04969
  * [24] Steenwinckel, B., Paepe, D.D., Hautte, S.V., Heyvaert, P., Bentefrit, M., Moens, P., Dimou, A., Bossche, B.V.D., Turck, F.D., Hoecke, S.V., Ongenae, F.: FLAGS: A methodology for adaptive anomaly detection and root cause analysis on sensor data streams by fusing expert knowledge with machine learning. Future Generation Computer Systems 116, 30–48 (Mar 2021). https://doi.org/10.1016/j.future.2020.10.015
  * [25] Tommasini, R., Della Valle, E., Mauri, A., Brambilla, M.: Rsplab: Rdf stream processing benchmarking made easy. In: d’Amato, C., Fernandez, M., Tamma, V., Lecue, F., Cudré-Mauroux, P., Sequeda, J., Lange, C., Heflin, J. (eds.) The Semantic Web – ISWC 2017. pp. 202–209. Springer International Publishing, Cham (2017). https://doi.org/10.1007/978-3-319-68204-4_21
  * [26] Zaharia, M., Xin, R.S., Wendell, P., Das, T., Armbrust, M., Dave, A., Meng, X., Rosen, J., Venkataraman, S., Franklin, M.J., Ghodsi, A., Gonzalez, J., Shenker, S., Stoica, I.: Apache spark: A unified engine for big data processing. Commun. ACM 59(11), 56–65 (oct 2016). https://doi.org/10.1145/2934664
  * [27] Zhang, Q., Song, Y., Routray, R.R., Shi, W.: Adaptive block and batch sizing for batched stream processing system. In: 2016 IEEE International Conference on Autonomic Computing (ICAC). pp. 35–44 (2016). https://doi.org/10.1109/ICAC.2016.27