Discretize-Then-Project (DtP)

• DtP is a design principle where discretization precedes projection, ensuring that critical discrete properties, such as pressure–velocity coupling, are preserved in reduced models.
• It is widely applied in areas like incompressible flow ROMs, dynamical low‐rank integrators, and quantum bath discretization, where projecting after discretization retains essential stability and conservation features.
• While DtP maintains structure and accuracy, it poses challenges including intrusiveness, operator scaling issues, and the need for tailored closures or hyper-reduction for complex modeling scenarios.

Searching arXiv for the cited DtP papers and closely related formulations. Discretize-Then-Project (DtP) denotes a class of constructions in which discretization precedes projection, reduction, or mapping. In reduced-order modeling for incompressible flow, DtP means first discretizing the full-order model and then projecting the resulting algebraic operators and vectors onto reduced spaces; in dynamical low-rank evolution, it means first discretizing the PDE in space to obtain a semi-discrete matrix ODE and then applying the projector-splitting integrator to that system; in quantum bath problems, it is a two-stage pipeline in which a continuum bath is first discretized into effective modes and then unitarily mapped to a finite chain (Star et al., 2020, Zhang et al., 21 Jul 2025, Vega et al., 2015, Rooholamin et al., 23 Jan 2026). Across these settings, DtP is used to preserve properties that are already encoded in the discrete model, rather than reconstructing them after a continuous projection.

1. Definition and relation to project-then-discretize

In the finite-volume reduced-order modeling literature on collocated grids, DtP is defined by building the reduced-order model after space discretization, or after both space and time discretization, so that the reduced operators are projections of the discrete full-order model rather than of the continuous PDE. For incompressible flow on collocated finite-volume grids, this distinction is consequential because the discrete pressure–velocity coupling depends on the precise algebraic pairing of divergence, gradient, and flux operators. The DtP formulation projects the actual discrete operators that already encode this coupling structure and the associated conservation properties (Star et al., 2020, Rooholamin et al., 23 Jan 2026).

The contrast with project-then-discretize (PtD) is explicit in the collocated-grid flow papers. PtD starts from the continuous PDE, projects it onto reduced spaces, and only then discretizes the reduced equations. In collocated finite-volume settings, PtD often fails to preserve the exact compatibility between discrete gradient and divergence and the specific stencil couplings needed to avoid checkerboard modes. The reported consequence is pressure–velocity decoupling after discretization unless auxiliary stabilization is added, for example supremizers, pressure penalty, or artificial compressibility (Rooholamin et al., 23 Jan 2026). By contrast, DtP avoids these additions by projecting the assembled discrete operators themselves.

In the dynamical low-rank setting studied for projector-splitting integrators (PSI), DtP likewise means discretizing once and then reusing the same discrete operator in all PSI substeps. PtD instead first projects onto the low-rank manifold in continuous space and then discretizes the projected PDEs, possibly with different spatial discretizations in different substeps. For the linear model problems analyzed in that work, DtP and PtD can nevertheless have the same stability restriction in some regimes, which shows that the DtP/PtD distinction is not uniformly equivalent to a stability advantage; its effect depends on the operator structure and the target problem class (Zhang et al., 21 Jul 2025).

2. DtP reduced-order models on collocated finite-volume grids

For incompressible Navier–Stokes equations on collocated finite-volume grids, DtP was developed as a Galerkin projection of the fully discrete full-order model. The discrete unknowns are cell-centered velocity $\boldsymbol{u}_p \in \mathbb{R}^{d h}$, cell-centered pressure $\boldsymbol{p}_p \in \mathbb{R}^h$, and face-centered velocity $\boldsymbol{u}_f \in \mathbb{R}^{d m}$. The finite-volume formulation uses the discrete continuity equation

$\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$

the interpolation

$$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$

and the semi-discrete momentum equation

$\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$

The explicit time integration is based on a Forward Euler projection method. Two variants are distinguished: the inconsistent flux method (IFM) and the consistent flux method (CFM) (Star et al., 2020).

The IFM computes pressure with a discrete pressure Poisson equation and updates the cell-centered velocity, but the face fluxes are not pressure-corrected. The CFM adds a face-velocity correction using a face-gradient operator $\boldsymbol{G}_f$,

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$

and, by construction, yields $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$. The paper emphasizes that this property is inherited at ROM level as well.

The reduced variables are introduced through

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$

with weighted orthonormality based on cell volumes and face areas. Projection of the fully discrete equations gives reduced pressure, momentum, and, for CFM, face-velocity equations. Because DtP projects the algebraic pressure equation and the boundary vectors, the ROM requires no pressure stabilization technique and no boundary control technique at ROM level. Pressure is computed from the same discrete pressure Poisson equation used by the full-order model, and boundary conditions are inherited by projection of the boundary vectors $\boldsymbol{p}_p \in \mathbb{R}^h$0, $\boldsymbol{p}_p \in \mathbb{R}^h$1, and $\boldsymbol{p}_p \in \mathbb{R}^h$2 (Star et al., 2020).

The reported numerical tests are a 2D lid-driven cavity and a 2D open cavity with inlet and outlet. In the lid-driven cavity, the domain is a square of length $\boldsymbol{p}_p \in \mathbb{R}^h$3 m on a structured $\boldsymbol{p}_p \in \mathbb{R}^h$4 grid, with $\boldsymbol{p}_p \in \mathbb{R}^h$5 m$\boldsymbol{p}_p \in \mathbb{R}^h$6/s, $\boldsymbol{p}_p \in \mathbb{R}^h$7, $\boldsymbol{p}_p \in \mathbb{R}^h$8 s, and $\boldsymbol{p}_p \in \mathbb{R}^h$9 s; with $\boldsymbol{u}_f \in \mathbb{R}^{d m}$0, velocity and pressure errors are approximately $\boldsymbol{u}_f \in \mathbb{R}^{d m}$1. In the open cavity, the grid has 7125 quadrilateral cells, $\boldsymbol{u}_f \in \mathbb{R}^{d m}$2 m$\boldsymbol{u}_f \in \mathbb{R}^{d m}$3/s, $\boldsymbol{u}_f \in \mathbb{R}^{d m}$4, $\boldsymbol{u}_f \in \mathbb{R}^{d m}$5 s, and $\boldsymbol{u}_f \in \mathbb{R}^{d m}$6 s; with $\boldsymbol{u}_f \in \mathbb{R}^{d m}$7, errors are approximately $\boldsymbol{u}_f \in \mathbb{R}^{d m}$8. The divergence behavior separates the two flux formulations: IFM yields local continuity errors of order $\boldsymbol{u}_f \in \mathbb{R}^{d m}$9 in the lid-driven cavity and $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$0 in the open cavity, whereas CFM yields errors of order $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$1. The open-cavity case also exhibits the largest reported speedups, with IFM ROM speedups of approximately $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$2 to $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$3 and CFM speedups of approximately $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$4 to $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$5 (Star et al., 2020).

3. Consistent flux DtP and hybrid turbulent-viscosity closure

The 2026 extension to turbulent incompressible flow on collocated grids preserves the DtP and consistent-flux construction but modifies the closure strategy for the turbulent viscosity field. The governing equations are

$\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$6

with LES Smagorinsky viscosity

$\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$7

and $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$8 for scalar transport in LES. The discrete formulation uses the same consistent-flux ingredients as the earlier collocated-grid work: face-based pressure gradients, the Laplacian $\boldsymbol{M}\,\boldsymbol{u}_f=\boldsymbol{0},$9, and a pressure Poisson equation that enforces discrete incompressibility at the face-flux level (Rooholamin et al., 23 Jan 2026).

The paper writes the reduced ansatz as

$$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$0

with reduced operators

$$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$1

and analogous projected convection, face-velocity, and source operators. The reduced pressure equation, reduced momentum equation, and reduced face-velocity equation retain the consistent-flux structure. Because the DtP ROM projects $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$2, $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$3, $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$4, $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$5, and $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$6 as assembled in the full-order model, the reported result is preservation of discrete compatibility and pressure–velocity coupling without Rhie–Chow interpolation at ROM level, pressure stabilization, or supremizers.

The distinctive feature of the turbulent formulation is the hybrid closure for $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$7. Direct intrusive projection of the eddy-viscosity field was found to be inconsistent in 3D convection-dominated regimes, so the paper adopts intrusive DtP for velocity and pressure but a non-intrusive data-driven closure for turbulent viscosity. The viscosity is represented as

$$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$8

and the temporal map is learned autoregressively with lookback window $$\boldsymbol{u}_f = \boldsymbol{I}_{p\rightarrow f}\,\boldsymbol{u}_p + \boldsymbol{u}_b,$$9,

$\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$0

Three neural architectures are evaluated: MLP, Transformer, and LSTM. The Transformer uses 4 heads, feed-forward dim 128, and dropout 0.1; the LSTM uses 64 then 32 units with return-sequences in the first layer. The training set comes from a 3D LES of a lid-driven cavity with 2000 snapshots, split into 1800 training and 200 validation samples; StandardScaler is used for inputs and outputs; the loss is mean squared error; Adam is the optimizer; the learning rates are $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$1 for LSTM, $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$2 for MLP, and $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$3 for Transformer; LSTM and MLP are trained for 1200 epochs, Transformer for 1000, with batch size 64 (Rooholamin et al., 23 Jan 2026).

The numerical experiment is a 3D lid-driven cavity in a unit cube on a structured $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$4 hexahedral mesh, with top-lid velocity $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$5, no-slip elsewhere, pressure zero-gradient on walls, and $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$6 at walls. The spatial discretization uses Gauss linear for gradients, linearUpwind for divergence, and linear orthogonal for Laplacian. Ten POD modes are retained for velocity, pressure, and turbulent viscosity. The reported best architecture is the LSTM, with typical relative errors of approximately $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$7 for velocity, approximately $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$8 for turbulent viscosity, approximately $\frac{\mathrm{d}\boldsymbol{u}_p}{\mathrm{d}t} = -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f)\,\boldsymbol{u}_p + \nu \boldsymbol{D}_p \boldsymbol{u}_p - \boldsymbol{G}_p \boldsymbol{p}_p + \boldsymbol{r}_p.$9 for pressure, approximately $\boldsymbol{G}_f$0 for energy, and approximately $\boldsymbol{G}_f$1 for enstrophy. The paper attributes the stability and low energy/enstrophy errors to the consistent-flux discretization inherited by the DtP ROM and the temporal memory of the LSTM closure (Rooholamin et al., 23 Jan 2026).

4. DtP in dynamical low-rank projector-splitting integrators

In the PSI literature, DtP is formulated for dynamical low-rank evolution after spatial discretization in one variable, so that the solution is collected in a matrix $\boldsymbol{G}_f$2 and approximated on the rank-$\boldsymbol{G}_f$3 manifold by

$\boldsymbol{G}_f$4

with $\boldsymbol{G}_f$5 and $\boldsymbol{G}_f$6. Given a matrix ODE

$\boldsymbol{G}_f$7

the projector-splitting integrator advances the solution through K-, S-, and L-steps. In the compact $\boldsymbol{G}_f$8-form used for DtP,

$\boldsymbol{G}_f$9

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$0

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$1

The defining DtP step is that the PDE is first discretized in space to obtain a semi-discrete matrix ODE $\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$2, and PSI is then applied directly to this system (Zhang et al., 21 Jul 2025).

For the linear hyperbolic model

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$3

the paper studies an upwind finite-difference discretization in $\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$4 and a Lie–Trotter PSI with Forward Euler substeps. The DtP amplification factor for a Fourier/eigen mode is

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$5

with

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$6

and the stability restriction

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$7

For the corresponding Lie–Trotter PtD formulation, the analysis gives the same CFL-type restriction. When Strang splitting with SSP-RK2 is used in the DtP substeps, the stability region is substantially enlarged; a sufficient condition found in the paper is

$\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$8

compared with $\boldsymbol{u}_f^{n+1} = \boldsymbol{I}_{p\rightarrow f}\boldsymbol{u}_p^n + \Delta t\,\boldsymbol{I}_{p\rightarrow f} \left( -\tilde{\boldsymbol{C}_p(\boldsymbol{u}_f^n)\boldsymbol{u}_p^n + \nu \boldsymbol{D}_p \boldsymbol{u}_p^n + \boldsymbol{r}_p \right) - \Delta t\,\boldsymbol{G}_f \boldsymbol{p}_p^{n+1},$9 for Lie–Trotter plus Forward Euler (Zhang et al., 21 Jul 2025).

For the linear parabolic model

$\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$0

the paper studies central differences in $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$1 and $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$2-schemes in time within the PSI substeps. If

$\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$3

the DtP amplification factor is

$\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$4

For $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$5 (Crank–Nicolson), the result is unconditional stability: $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$6 For $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$7 (Backward Euler), a sufficient condition is

$\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$8

For $\boldsymbol{M}\boldsymbol{u}_f^{n+1}=0$9 (Forward Euler), a sufficient condition is

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$0

The paper also analyzes the negative S-step phenomenon: the S-step corresponds to marching backward in diffusion time, but with Crank–Nicolson in each substep the overall amplification collapses to the stable Crank–Nicolson factor. A hybrid scheme using Backward Euler in K- and L-steps and Forward Euler in the S-step yields

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$1

and is therefore unconditionally stable. For the linear diffusion operator with central differences, DtP and PtD have identical stability properties (Zhang et al., 21 Jul 2025).

5. DtP for quantum bath discretization and chain mapping

In quantum bath problems, DtP refers to a different projection target. The continuous system-plus-bath Hamiltonian is written in star form as

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$2

with

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$3

A discrete approximation with $$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$4 bath modes is

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$5

The spectral density is

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$6

and the hybridization function is

$$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$7

The DtP workflow is then: first discretize the continuum bath into effective modes $$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$8, and then project or map the discrete star Hamiltonian onto a finite one-dimensional tight-binding chain suitable for real-time methods such as TEBD or DMRG (Vega et al., 2015).

The paper classifies discretization strategies into direct discretization, orthogonal polynomial (OP), and numerical optimization. Direct discretization chooses intervals $$\boldsymbol{u}_p \approx \boldsymbol{\Phi}\,\boldsymbol{a},\qquad \boldsymbol{p}_p \approx \boldsymbol{X}\,\boldsymbol{b},\qquad \boldsymbol{u}_f \approx \boldsymbol{\Psi}\,\boldsymbol{c},$$9 and defines

$\boldsymbol{p}_p \in \mathbb{R}^h$00

The OP strategy uses Gauss–Christoffel quadrature with weight $\boldsymbol{p}_p \in \mathbb{R}^h$01, either $\boldsymbol{p}_p \in \mathbb{R}^h$02 or $\boldsymbol{p}_p \in \mathbb{R}^h$03. The orthonormal polynomials satisfy

$\boldsymbol{p}_p \in \mathbb{R}^h$04

and the recurrence

$\boldsymbol{p}_p \in \mathbb{R}^h$05

For $\boldsymbol{p}_p \in \mathbb{R}^h$06, the resulting chain Hamiltonian is

$\boldsymbol{p}_p \in \mathbb{R}^h$07

with

$\boldsymbol{p}_p \in \mathbb{R}^h$08

and coupling

$\boldsymbol{p}_p \in \mathbb{R}^h$09

The OP mapping with $\boldsymbol{p}_p \in \mathbb{R}^h$10 is equivalent to chain mapping via Lanczos tridiagonalization (Vega et al., 2015).

For quadratic Hamiltonians, the paper states that the OP strategy is the best strategy in the sense that it gives the numerically exact time evolution up to a maximum time $\boldsymbol{p}_p \in \mathbb{R}^h$11. The formulas are

$\boldsymbol{p}_p \in \mathbb{R}^h$12

where $\boldsymbol{p}_p \in \mathbb{R}^h$13 is the bath bandwidth. The interpretation given is finite-size revival in the finite chain: before $\boldsymbol{p}_p \in \mathbb{R}^h$14, Gauss–Christoffel quadrature integrates the relevant polynomial expansion exactly up to degree $\boldsymbol{p}_p \in \mathbb{R}^h$15. For non-quadratic Hamiltonians, the paper states that no such best strategy exists. In that regime, direct discretizations that place nodes where $\boldsymbol{p}_p \in \mathbb{R}^h$16 is large can outperform BSDO/OP for target observables, as illustrated by a DMFT-motivated SIAM bath with three Gaussian peaks and gaps. The paper also reports that numerical optimization on the imaginary axis is stable only for small $\boldsymbol{p}_p \in \mathbb{R}^h$17, approximately up to 15, and is not a controlled real-time strategy (Vega et al., 2015).

6. Common properties, advantages, and limitations

Across the supplied literature, DtP is consistently used to preserve structure that is already present in the discrete representation. In collocated-grid incompressible-flow ROMs, this means preserving the assembled pressure–velocity coupling, the boundary treatment, and, in the consistent-flux setting, mass conservation at the face-flux level without ROM-level pressure stabilization, supremizers, pressure penalty, artificial compressibility, or boundary control (Star et al., 2020, Rooholamin et al., 23 Jan 2026). In PSI, it means that stability can be analyzed directly from the discrete operator symbol, yielding closed-form bounds for hyperbolic and parabolic model problems and clarifying the role of the negative S-step (Zhang et al., 21 Jul 2025). In quantum bath discretization, it means that exactness up to a computable time horizon can be tied to the finite chain obtained after discretizing the continuum bath and mapping it with orthogonal polynomials (Vega et al., 2015).

The limitations are equally domain-specific. In the collocated finite-volume ROM setting, DtP is intrusive, depends on access to discrete operators and boundary vectors, and inherits the discretization choices of the full-order model; convection tensors scale cubically with reduced dimension, and hyper-reduction such as DEIM or GNAT is proposed as an extension (Star et al., 2020, Rooholamin et al., 23 Jan 2026). In the turbulent extension, the neural turbulent-viscosity closure must be trained on representative data, may not extrapolate well to out-of-distribution regimes or parameter changes without retraining, and may require physical constraints or hyper-reduction when online operator dependence on $\boldsymbol{p}_p \in \mathbb{R}^h$18 is strong (Rooholamin et al., 23 Jan 2026). In PSI, Lie–Trotter plus Forward Euler is restrictive for hyperbolic problems, while stiff diffusion favors Crank–Nicolson or the hybrid BE–FE scheme; the linear analysis also assumes periodic boundary conditions, constant coefficients, and exact orthonormality (Zhang et al., 21 Jul 2025). In the bath-discretization setting, the optimality guarantee is limited to quadratic Hamiltonians; for interacting systems, observable convergence rather than hybridization fidelity alone becomes decisive, and no universal best discretization exists (Vega et al., 2015).

A plausible implication is that DtP is less a single algorithm than a design principle: whenever critical conservation, coupling, or stability properties are encoded by the discrete model rather than the continuous equation alone, projecting after discretization can reproduce those properties more faithfully than projecting before discretization. The supplied works support that interpretation in three distinct senses: exact replication of discrete pressure–velocity coupling on collocated grids, stability characterization of low-rank splitting on discrete operators, and controlled real-time fidelity of discrete chain representations for quadratic baths (Star et al., 2020, Zhang et al., 21 Jul 2025, Vega et al., 2015, Rooholamin et al., 23 Jan 2026).