“…In essence, it is to explore whether or not the range of the high-dimensional function E h is always contained in G: E h (G s+k+1 ) ⊆ G. For some scalar PDEs with linear constraints, for instance, the scalar conservation laws with the constraints linearly defined by maximum principle, a general approach for bound-preserving analysis and design is to exploit certain monotonicity in schemes; see, e.g., [67,14,31]. Yet, for PDE systems especially with nonlinear constraints, there is no unified tool like monotonicity, so that direct and complicated algebraic verification usually has to be performed for each constraint case-by-case for different schemes and different PDEs; see, e.g., [68,38,56,41,66,35,53]. Therefore, the design and analysis of bound-preserving schemes involving nonlinear constraints are highly nontrivial, even for first-order schemes; cf.…”