  We introduce a fully generalized quiescent chemical reactor system in
arbitrary space $\vdim =1,2$ or 3, with $n\in\mathbb{N}$ chemical constituents
$\alpha_{i}$, where the character of the numerical solution is strongly
determined by the relative scaling between the local reactivity of species
$\alpha_{i}$ and the local functional diffusivity $\mathscr{D}_{ij}(\alpha)$ of
the reaction mixture. We develop an operator time-splitting predictor
multi-corrector RK--LDG scheme, and utilize $hp$-adaptivity relying only on the
entropy $\mathscr{S}_{\mathfrak{R}}$ of the reactive system $\mathfrak{R}$.
This condition preserves these bounded nonlinear entropy functionals as a
necessarily enforced stability condition on the coupled system. We apply this
scheme to a number of application problems in chemical kinetics; including a
difficult classical problem arising in nonequilibrium thermodynamics known as
the Belousov-Zhabotinskii reaction where we utilize a concentration-dependent
diffusivity tensor $\mathscr{D}_{ij}(\alpha)$, in addition to solving a simple
equilibrium problem in order to evaluate the numerical error behavior.
