## 3 Numerical Framework
### Implicit Time Integration
Notice that all the terms appearing in (21) are well defined, in particular, regardless of the parameter \(\gamma\), \(\left(\mathbf{M}_{P}\right)_{\mathrm{k}}^{\mathrm{i}}(\mathbf{M}_{P}^{-1})^{ \mathrm{j}}\Big{|}_{\mathrm{K}}=\)  
\[=\begin{bmatrix}\mathbf{I}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\left(\mathbf{B}_{\rho}\right)_{\mathrm{k}}^{\mathrm{i}}\left( \mathbf{B}_{\rho}^{-1}\right)^{\mathrm{j}}&\mathbf{0}&-\left(\mathbf{B}_{\rho} \right)_{\mathrm{k}}^{\mathrm{i}}\left(\mathbf{B}_{\rho}^{-1}\mathbf{B}_{\rho \bar{\rho}^{\mathrm{i}}\mu_{1}}\mathbf{B}_{\rho\bar{\rho}}^{-1}\right)^{ \mathrm{j}}+\left(\mathbf{B}_{\rho\bar{\rho}\mu_{1}}\right)_{\mathrm{k}}^{ \mathrm{i}}\left(\mathbf{B}_{\bar{\rho}}^{-1}\right)^{\mathrm{j}}\\ \mathbf{0}&\mathbf{0}&\left(\mathbf{B}_{\rho}\right)_{\mathrm{k}}^{\mathrm{i}} \left(\mathbf{B}_{\rho}^{-1}\right)^{\mathrm{j}}&-\left(\mathbf{B}_{\rho} \right)_{\mathrm{k}}^{\mathrm{i}}\left(\mathbf{B}_{\rho}^{-1}\mathbf{B}_{\rho \bar{\rho}^{\mathrm{i}}\mu_{2}}\mathbf{B}_{\rho\bar{\rho}^{-1}}^{-1}\right)^{ \mathrm{j}}+\left(\mathbf{B}_{\rho\bar{\rho}\mu_{2}}\right)_{\mathrm{k}}^{ \mathrm{i}}\left(\mathbf{B}_{\bar{\rho}}  
minimizing the computational effort to achieve a user-defined accuracy and (ii) improve the overall time integration robustness by limiting the LTE.  
Since the LTE cannot be computed directly, as the analytical solution of the problem is (in general) unknown, an LTE estimation procedure is required. This task is usually carried out by comparing the numerical solution \(\mathbf{W}^{n+1}\) with a less-accurate one \(\hat{\mathbf{W}}^{n+1}\). Accordingly, the LTE estimator is defined as  
\[r^{n+1}=\|\mathbf{W}^{n+1}-\hat{\mathbf{W}}^{n+1}\|,\]  
where \(\|\cdot\|\) is a user-defined norm. At each time step of an ESDIRK scheme, the less-accurate/lower order solution, commonly referred to as the _embedded_ solution, can be computed by replacing the set of coefficients \(b\) in Equation (18), with a modified set \(\hat{b}\), as follows  
\[\hat{\mathbf{W}}^{n+1}=\mathbf{W}^{n}-\Delta t^{n}\sum_{i=1}^{\mathbf{s}}\hat{ b}_{i}\left(\mathbf{M}_{P}^{-1}\mathbf{R}\right)^{i}.\]  
At the end of each time step, the time step adaptation algorithm of [37,38] checks if the following condition is satisfied  
\[r^{n+1}<\zeta\eta, \tag{22}\]  
where the threshold value \(\eta\) of the LTE estimator and the safety factor \(\zeta\) are user-defined parameters. If (22) is verified, then the solution \(\mathbf{W}^{n+1}\) is accepted and the size of the next time step is set as  
\[\Delta t^{n+1}=\left(\frac{\eta}{\psi^{n+2}}\right)^{\frac{1}{q}}; \tag{23}\]  
otherwise, \(\mathbf{W}^{n+1}\) is rejected and the time step is recomputed, imposing  
\[\Delta t^{n}=\left(\frac{\eta}{\psi^{n+1}}\right)^{\frac{1}{q}}. \tag{24}\]  
The parameter \(q\) is the order of accuracy of the LTE estimator, which depends on the particular combination of ESDIRK scheme/embedded solution considered. The estimator constant \(\psi\) at the end of each step can be easily computed by exploiting the asymptotic behavior  
\[\psi^{n+1}=\frac{r^{n+1}}{\left(\Delta t^{n}\right)^{q}}. \tag{25}\]  
Notice that in Equation (24) the constant \(\psi^{n+1}\) is computed in Equation (25) using the size of the rejected time step. Conversely, the future time value \(\psi^{n+2}\) required in Equation (23) is unknown and must be extrapolated using information from previous time levels. In this work, the extrapolation is based on the H211b controller, see [39,40] for further details. In order to improve the robustness of the adaptive algorithm, a 'bound to the maximum time step size increase/reduction' is enforced by a smooth limiter function,  
\[\overline{\Delta t}^{*}=\Delta t^{n}\bigg{[}1+\zeta\arctan\left(\frac{ \Delta t^{*}-\Delta t^{n}}{\xi\Delta t^{n}}\right)\bigg{]},\]  
where \(\Delta t^{*}\) is computed by either Equation (23) or Equation (24) in agreement with the condition (22) and \(\overline{\Delta t}^{*}\) is the new limited value for the time step size.