## 1 Introduction  
Adaptive numerical methods for solving ordinary differential equations (ODEs) have become a significant topic in numerical analysis. The primary motivation for this approach lies in the challenges posed by stiff ODEs, where sudden changes in certain components may occur, and where explicit methods are subject to severe stability restrictions. Such behavior is common, for example, in models of chemical reactions that involve higher-order reaction terms. The steep gradients necessitate very fine time-stepping, which can substantially slow down the simulation. It is important to note that realistic chemical reaction systems often involve a large number of reaction steps with widely varying reaction-rate coefficients [1]. Similarly, the semi-discretization of large-scale air pollution models [2] can result in ODE systems with a vast number of unknowns. To address this, it is essential to identify regions where steep gradients occur, applying small time steps in their vicinity. Several algorithms have been proposed in the literature to construct efficient adaptive refinements (see [3], Section 5.2). We begin by noting that any spatial semi-discretization of partial differential equations (PDEs) yields a large, stiff system of ordinary differential equations (ODEs). As a result, efficient numerical methods are required--even when the reaction terms are nonstiff. However, heuristic algorithms for selecting time steps can be prone to instabilities. For instance, it is important to adjust the time steps smoothly by controlling the ratio between successive step sizes. This necessity naturally connects time-step selection with control theory. Indeed, some of the most effective adaptive time-stepping algorithms have been developed using control-theoretic approaches. Further details can be found in the classical works [4; 5; 6; 7]. For a recent approach, we refer readers to [8]. Another strategy for accelerating the numerical solution of ODEs is to use parallel-in-time algorithms. Over the past decade, interest in these methods has been revitalized by the advent of massively parallel computing architectures. In addition to traditional applications involving PDEs [9; 10; 11; 12], recent developments such as neural ordinary differential equations [13] and physics-informed neural networks [14] have further broadened the scope of parallel-in-time methods, especially in machine learning and deep learning contexts [15; 16; 17; 18]. Parallel-in-time methods can be broadly categorized into four main groups: multiple shooting methods; domain decomposition and waveform relaxation methods; space-time multigrid methods; and direct parallel-in-time algorithms. A comprehensive overview of these techniques is provided in the review paper [19] and the monograph [20]. In this work, we focus on the multigrid reduction-in-time (MGRIT) method introduced in [21]. MGRIT is closely related to the parareal algorithm developed in [22] and has since been extended in a variety of directions. These include applications to nonlinear parabolic equations [23], adaptive (variable-step) multistep methods [24], layer-parallel training of deep residual neural networks [25], efficiency improvements via Richardson extrapolation [26], and hyperbolic problems such as linear advection [27]. Adaptive parallel-in-time methods have been explored in [24; 26], with a primary focus on metrics such as time versus discretization error and processor versus wall-clock time. However, another important aspect of reliability is how the global error responds to changes in the specified tolerance. Ideally, for problems with regular solutions, a small change in the tolerance should result in a proportionally small change in the error. This property, referred to as tolerance proportionality, is fundamental to computational stability [28].