{"id":"6d3e0537-4d6a-4b37-a8b7-810c9e731e85","arxiv_id":"2411.14315","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"A stabilized finite element harmonic balance solver for the Navier-Stokes equations achieves one to two orders of magnitude speedup for periodic cardiovascular flow simulations with good accuracy.","lead":"This paper presents a harmonic balance (frequency-domain) finite element solver for periodic blood flow simulations that computes all time points in a cardiac cycle simultaneously. The solver is reported to run 30 to 160 times faster than a conventional time-stepping solver while matching its results to within a few percent on three patient-specific cases.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (7) is not a direct transform of the anti-aliased frequency equations: the pointwise convective term reintroduces the aliasing that A_i's zero corners were meant to exclude.","rationale":"The reader correctly identified the false diagonalization E^{-1} A_i E = diag(u_i) as the weakest assumption, and this is indeed the central theoretical flaw. My analysis extends that concern: the false diagonalization is not just a cosmetic derivation error, because the two formulations have different aliasing behavior. The frequency-domain Eq. (4) uses a truncated convolution that discards product modes above M-1, while the time-domain Eq. (7) uses a pointwise product at N = 2M-1 collocation points, which aliases those high-frequency products into the resolved band. Therefore, Eq. (7) cannot be a direct transformation of Eq. (4), and the paper's claim that aliasing is eliminated is internally inconsistent with the stated U_i = diag(u_i). This is a genuine correctness risk for the method's theoretical foundation. However, the paper's empirical validation against the time-stepping solver provides independent support for the practical accuracy of the implemented collocation formulation on the three test cases. The reported speedups are also plausible and are backed by wall-clock measurements. For these reasons, the conditional verdict remains appropriate: the method may work, but the paper's derivation and its anti-aliasing claim need correction or clarification, and a numerical comparison between the aliased and anti-aliased convective treatments would settle whether the discrepancy actually affects accuracy in practice.","tokens_in":24018,"tokens_out":14984,"duration_ms":143038,"concrete_test":"For a small case (e.g., the N = 7 pulmonary or cerebral simulation), compute the convective term at each time point using two routes: (i) the implemented pointwise product U_j ∂u_i/∂x_j, and (ii) the frequency-domain product A_j ∂u_i/∂x_j (with the zero-corner truncation) transformed back to time via E. If the two residuals differ by more than the reported convergence tolerance, Eq. (7) is not a transformation of Eq. (4). Then re-run the N = 7 case with the anti-aliased convective term (option ii) and compare the velocity error versus the time solver; if the error drops significantly relative to the implemented version, the aliasing error is load-bearing for the accuracy claim.","verdict_should_be":"UNCHANGED","load_bearing_attack":"In Section 2.1, the paper defines A_i in Eq. (4) with zero entries at the top-right and bottom-left corners precisely to 'eliminate the aliasing effect' by preventing product frequencies at or above M from entering the resolved band. The derivation of Eq. (7) claims E^{-1} A_i E = diag(u_i), which would require A_i to be circulant. A_i is not circulant because of those zero corners, so the diagonalization is false. The consequence is not merely a missing proof: the convective term ρ U_j ∂u_i/∂x_j with U_j = diag(u_i) is a pointwise product sampled at N = 2M-1 points. The discrete product of two M-mode fields contains frequencies up to 2M-2, and with N = 2M-1 samples those high frequencies alias into the resolved modes when the collocation residual is formed. Thus the harmonic balance equations (7) do not inherit the anti-aliasing property of Eq. (4); they correspond to a different, aliased convolution. Remark 1's alternative derivation from Eq. (1) yields a valid collocation method, but it does not establish equivalence with Eq. (4) or justify the claim that aliasing is eliminated. The reported accuracy for smooth, low-to-moderate Reynolds flows may be insensitive to this aliasing because the high-mode energy is small, but the paper's theoretical foundation for the method is incorrect as stated.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents a stabilized finite element solver for the incompressible Navier-Stokes equations in harmonic balance form, targeting time-periodic cardiovascular flows. The unknowns at N uniformly spaced instants in the cardiac cycle are solved simultaneously; the time-derivative coupling is implemented through an FFT, and matrix splitting is used to make the cost scale roughly linearly with N. The solver is validated on three patient-specific geometries (Glenn pulmonary flow, cerebral arteries, left main coronary artery) by comparing velocity, pressure, and outlet flow rates against a conventional GLS time-stepping solver on the same meshes. Reported speedups range from about 10 to 100 times, with integrated velocity errors below 5% when N is sufficiently large (N >= 19 in the cases shown), and memory overhead of about 2.4x at N = 25. The authors conclude that the solver is viable for low-to-intermediate Reynolds number periodic flows.","tokens_in":24287,"tokens_out":11582,"duration_ms":114676,"significance":"The numerical contribution is potentially valuable: it addresses a practical bottleneck in patient-specific cardiovascular simulation, demonstrates near-linear cost scaling through FFT-based time coupling, and reports careful comparisons on three clinically relevant geometries using the same meshes, stabilization framework, and hardware for both solvers. The accuracy assessment is benchmarked against an independent time-stepping solver rather than being self-referential, which is a genuine strength. However, the theoretical derivation connecting the frequency-domain formulation to the harmonic balance form contains an invalid algebraic step, and the manuscript overstates the equivalence between Eq. (4) and Eq. (7). The empirical results may survive a corrected derivation, but the formulation's mathematical foundation and its anti-aliasing claims need to be revised before the paper is fully sound.","major_comments":[{"comment":"The derivation of Eq. (7) from Eq. (4) implicitly asserts that E^{-1} A_i E = diag(u_i) for the convective coupling matrix. This assertion is false. The matrix A_i defined after Eq. (4) is a truncated Toeplitz matrix with zero entries at its top-right and bottom-left corners, enforced by the condition |j-k| < M to 'eliminate aliasing'; it is not a circulant matrix, and the DFT diagonalizes only circulant matrices. The pointwise product U_j \\partial u_i / \\partial x_j in Eq. (7) corresponds to a circular convolution that wraps contributions from modes |m| >= M back into the resolved band, so Eq. (7) does not inherit the anti-aliasing property claimed for Eq. (4). Remark 1 correctly observes that Eq. (7) can be derived from Eq. (1) as a spectral collocation method, but that alternative derivation does not establish equivalence with Eq. (4) and does not justify the claim that aliasing is eliminated. Because the paper's central methodological premise is that Eq. (7) is a direct transform of the anti-aliased frequency equations, this step must be corrected or the method must be explicitly reframed as a collocation formulation.","section":"Section 2.1, Eqs. (4)-(7)"},{"comment":"There is an indexing inconsistency in the transform definitions. The vector u*_i in Eq. (3) is ordered by mode index -M+1, ..., M-1, whereas the DFT matrix defined in Eq. (5), E(j,k) = (1/N)e^{-2\\pi i j k/N} with j,k = 0, ..., N-1, orders its rows by DFT frequency 0, ..., N-1. As written, the equality u*_i = E u_i does not hold for the stated ordering, and H = E^{-1} \\Omega E is not the time-derivative matrix for the stated unknowns unless an intermediate permutation is introduced. The implementation may use a standard FFT ordering internally, but the manuscript's formulation is not self-consistent without explicitly specifying this permutation.","section":"Section 2.1, Eqs. (3)-(5)"},{"comment":"The central speedup claims (10 to 100 times faster) are measured against a conventional time-stepping solver that was run for four cardiac cycles, but no cycle-to-cycle convergence metric is reported. If the time-stepping reference had not reached the time-periodic state to within the reported 5% error levels, the speedup could be inflated and the error integrals would not have a converged reference. Please add a periodic-convergence diagnostic for the time solver, such as the change in velocity or pressure norms between successive cardiac cycles, or state explicitly the criterion used to consider the reference solution converged.","section":"Section 3, speedup comparisons"}],"minor_comments":[{"comment":"The caption reads 'R V A, left vertebral artery'; it should read 'right vertebral artery'.","section":"Figure 8 caption"},{"comment":"The claim that the equivalent time step 2\\pi/(N\\omega) is 'large enough' to omit the (1/\\Delta t)^2 term from the stabilization parameter should be quantified; for N = 25 the equivalent time step is not obviously large compared with the convective time scale, and reporting the element Womersley number for each case would make the claim checkable.","section":"Section 2.2, Remark 4"},{"comment":"The terminology is inconsistent: the paper refers to the 'harmonic balanced method' and the 'harmonics-balance Navier-Stokes equations' in places; the standard term 'harmonic balance' should be used consistently.","section":"Abstract and Section 1"},{"comment":"The data availability statement says results are available on request; depositing the meshes, boundary conditions, and scripts in a permanent repository would improve reproducibility.","section":"Data availability"}],"recommendation":"major_revision","confidential_remarks":"The invalid diagonalization step in Section 2.1 is the main technical concern; it is fixable by rewriting the derivation as a time-spectral collocation method or by explicitly proving a corrected equivalence, and the independent numerical validation suggests the empirical conclusions may stand. The indexing inconsistency between the frequency-vector ordering and the DFT matrix should also be resolved in revision. The paper's relation to the authors' prior frequency-domain solver [31] and to the harmonic balance formulation of Hall et al. [47] should be clarified, since the novelty claim rests on the linear-cost scaling rather than on the harmonic balance idea itself."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Here's my take on the harmonic balance solver paper. The headline result is real: the GLS-stabilized harmonic balance solver runs 30 to 160 times faster than a conventional time-stepping solver on three patient-derived geometries, with velocity errors under 5 percent when the mode count is high enough (N = 19 or 25). The numerical validation is careful and honest. Both solvers use the same GLS formulation; the error is measured against the time solver, which is the right benchmark; and the paper openly discusses the cases it cannot handle, like high-Reynolds aorta flow.\n\nWhat's genuinely new is the first stabilized finite element implementation of harmonic balance Navier-Stokes with FFT-based matrix-vector products, giving roughly linear cost scaling in the number of time levels. The authors demonstrate convergence, strong scaling, and memory use. That is a useful step for cardiovascular simulation, and the performance claims are plausible from the reported experiments.\n\nThe soft spot is the derivation in Section 2.1. The paper states that Eq. (7) is a direct transformation of the frequency equations (4), via E^{-1} A_i E = diag(u_i). That is false. The matrix A_i has zero corners by construction to eliminate aliasing, so it is not circulant, and the DFT does not diagonalize it. Consequently Eq. (7) is not the transform of the anti-aliased equations. What Eq. (7) actually is, as Remark 1 suggests, is a spectral-collocation discretization of the original time-domain equations. That is a legitimate method and the numerical results probably survive the correction, because for smooth low-Reynolds flows the aliasing from the pointwise product is small. But the theoretical claim as written is wrong, and the paper's own emphasis on 'eliminating the aliasing effect' is not inherited by the scheme it actually solves. This should be fixed, and the authors should either prove or empirically show that the aliasing is harmless for these flows. On reproducibility, 'data available on request' is too weak for a computational paper making orders-of-magnitude speedup claims; code and input files would make the results checkable.\n\nBottom line: the paper is solid where it matters most, the empirical results appear trustworthy, and the derivation flaw is identifiable and fixable. This deserves a serious referee and a revision, not a desk rejection. I'd bring it to our reading group to discuss the aliasing question, and I'd cite it if I worked on time-spectral methods.","headline":"Useful practical contribution with a flawed derivation of the harmonic balance equations; the numerical evidence is solid and the flaw is fixable, so it deserves review.","tokens_in":24810,"tokens_out":3297,"would_cite":true,"duration_ms":30438,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65M60","65T50","76D05","76Z05"],"pacs":[],"model":"deepseek-v4-flash","headline":"A harmonic balance finite element solver claims to compute time-periodic cardiovascular flows 10 to 100 times faster than conventional time-stepping, with under 5 percent velocity error at 19 modes.","keywords":["harmonic balance","stabilized finite element method","cardiovascular simulation","incompressible Navier-Stokes equations","time-periodic flow","frequency domain","fast Fourier transform","patient-specific hemodynamics"],"falsifier":"Compute $E^{-1}A_iE$ for a nonzero velocity field with $M \\geq 3$ and inspect the off-diagonal entries; if any are nonzero, the claimed diagonalization behind the $O(N\\log N)$ cost is not exact. A stronger test is to run the harmonic balance solver with the full convection coupling matrix $A_i$, without the diagonalizing step, and compare residuals and converged solutions against the diagonal version on a case with $N \\geq 19$; equal results would show the omitted couplings are negligible, while different results would identify the error introduced by the diagonalization.","tokens_in":23767,"feed_emoji":"🫀","tokens_out":6837,"duration_ms":59822,"temperature":0.7,"pith_summary":"Cardiovascular blood flow is periodic, so the authors argue that a spectral discretization in time should replace thousands of time steps with a handful of Fourier modes. They build a stabilized finite element solver for the harmonic balance form of the incompressible Navier-Stokes equations, in which all time points are solved simultaneously and the time-derivative term is a coupling matrix computed by fast Fourier transforms. On three patient-specific geometries, a Glenn anastomosis, a cerebral artery, and a left main coronary artery, the solver reproduces the conventional time-stepping solution with velocity errors below about 5 percent once the number of modes is large enough, while cutting wall-clock time by one to two orders of magnitude. The practical stake is that patient-specific hemodynamic simulations, which currently take many hours to days, could become fast enough for clinical decision-making.","feed_headline":"Harmonic balance solver cuts blood-flow simulation time 10-100x","feed_subtitle":"Periodic cardiovascular flows solved in frequency space match time-stepping within about 5% velocity error at N=19.","key_machinery":"The central object is the harmonic balance form of the Navier-Stokes equations, obtained by writing the frequency-domain equations with unknowns organized as discrete Fourier coefficient vectors and applying the inverse discrete Fourier transform matrix $E^{-1}$. The load-bearing algebraic objects are the time-coupling matrix $H = E^{-1}\\Omega E$, a skew-symmetric zero-diagonal matrix coupling all $N$ time points that is applied via fast Fourier transforms, and the diagonal matrices $U_j = \\operatorname{diag}(u_j)$ for the convective term. The paper's cost scaling depends on the claim that the discrete Fourier transform diagonalizes the convective coupling, leaving a diagonal $U_j$; this is what allows the tangent matrix to split into a pointwise block $P$ plus an FFT-applied coupling block $C$, and it is what makes the residual and tangent computations approximately $O(N\\log N)$ rather than $O(N^2)$. The stabilized weak form uses Galerkin/least-squares with a stabilization parameter $\\tau$ chosen to recover the steady-state form and to keep the method consistent as the equivalent time step vanishes.","core_discovery":"The paper's central claim is that the harmonic balance form of the incompressible Navier-Stokes equations, discretized with a Galerkin/least-squares stabilized finite element method, is an accurate and much faster replacement for conventional time-stepping in physiologically stable, time-periodic blood flows. The authors represent velocity, pressure, and boundary data as $N$ equally spaced time points per cardiac cycle, obtained from a truncated discrete Fourier series, and solve one coupled system containing all time points. The time-coupling matrix $H = E^{-1}\\Omega E$, where $\\Omega$ contains the imaginary angular frequencies and $E$ is the discrete Fourier transform matrix, is multiplied using fast Fourier transforms, giving an approximately $O(N\\log N)$ cost instead of the $O(N^2)$ cost of their earlier frequency-domain solver. In the three test cases, the harmonic balance solver matched the time-stepping solver's velocity field within roughly 5 percent at $N = 19$ and produced speedups from 30-fold to 40-fold for the largest $N$, with the paper reporting a 10-to-100-fold speedup across all cases. The paper also finds that the truncation error of the inlet boundary condition predicts the overall solution error at low Reynolds number, and that the harmonic balance solution acts as a built-in filter for noisy clinical inflow data.","pith_inferences":["Editorial inference: the claimed diagonalization $E^{-1}A_iE = \\operatorname{diag}(u_i)$ in Section 2.1 is not generally true because $A_i$ has zero top-right and bottom-left corners and is not circulant; if off-diagonal couplings remain, the linear-cost scaling and the direct-transform claim would need to be revisited, although the numerical results could still hold if those couplings are small i","Editorial inference: the same harmonic balance machinery could apply to other periodic physiological flows, such as respiratory airflow, where the period is known and Reynolds numbers are comparable; the paper mentions this but does not test it.","Editorial inference: the built-in spectral truncation could be used deliberately as a denoising step for clinical inflow waveforms, trading a controlled truncation error for robustness against measurement noise."],"forward_implications":["For flows whose boundary conditions are well represented by 19 modes, the harmonic balance solver matches conventional time-stepping within about 5 percent velocity error, and within 3 percent for outlet flow rates even at $N = 7$ to $N = 13$.","The approximately linear cost scaling removes the quadratic scaling that limited the earlier frequency-domain approach, making large mode counts practical.","Speedups of 30 to 40 times at the largest $N$ tested mean a 30-hour coronary simulation can be completed in roughly 30 minutes with similar accuracy.","Memory overhead stays modest, about 2.4 times the time-stepping solver at $N = 25$, because most memory is occupied by sparse matrices and connectivity rather than unknowns.","Because the equations are presented as a direct transform of the frequency equations, the method can be adapted into an existing stabilized finite element Navier-Stokes code structure with relatively small changes."],"supporting_citations":[{"why":"Supplies the harmonic balance form of the Navier-Stokes equations that the paper transforms and implements.","marker":"[47]"},{"why":"Prior frequency-domain stabilized finite element solver that this work extends, providing the GLS stabilization and pseudo-time-stepping scheme.","marker":"[31]"},{"why":"Introduces the augmented stabilization method and analysis used to justify the GLS approach in the time-spectral setting.","marker":"[30]"},{"why":"FFTW library implements the fast Fourier transforms that give the $O(N\\log N)$ matrix-vector products for the time-coupling matrix.","marker":"[65]"},{"why":"Supplies the patient-specific geometries used in the three test cases.","marker":"[76]"},{"why":"Provides the Glenn procedure boundary conditions and clinically relevant simulation setup.","marker":"[77]"},{"why":"Provides the coronary inflow profile and prior CFD setup for the left main coronary case.","marker":"[86]"}],"fun_headline_variants":["Harmonic balance solver cuts heart flow sim time 100x","Blood flow modeling 100x faster with harmonic balance","Frequency-space heart simulation runs in minutes, not hours","New harmonic balance method accelerates cardiovascular sims 10-100x","Heart flow simulations in 30 minutes via harmonic balance"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The derivation of the harmonic balance equations from the frequency-domain equations in Section 2.1 assumes that the discrete Fourier transform diagonalizes the convective coupling matrix $A_i$, i.e. $E^{-1}A_iE = \\operatorname{diag}(u_i)$, but $A_i$ is not circulant because of its zero top-right and bottom-left corners, so this diagonalization does not hold in general.","fun_headline_variants_meta":{"raw":{"variants":["Harmonic balance solver cuts heart flow sim time 100x","Blood flow modeling 100x faster with harmonic balance","Frequency-space heart simulation runs in minutes, not hours","New harmonic balance method accelerates cardiovascular sims 10-100x","Heart flow simulations in 30 minutes via harmonic balance"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000787,"raw_usage":{"total_tokens":3532,"prompt_tokens":1065,"completion_tokens":2467,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":681,"completion_tokens_details":{"reasoning_tokens":2386}},"tokens_in":681,"tokens_out":2467,"duration_ms":17006,"temperature":1.0,"reasoning_tokens":2386,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T15:19:42.107182+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute $E^{-1}A_iE$ for a nonzero velocity field with $M \\geq 3$ and inspect the off-diagonal entries; if any are nonzero, the claimed diagonalization behind the $O(N\\log N)$ cost is not exact. A stronger test is to run the harmonic balance solver with the full convection coupling matrix $A_i$, without the diagonalizing step, and compare residuals and converged solutions against the diagonal version on a case with $N \\geq 19$; equal results would show the omitted couplings are negligible, while different results would identify the error introduced by the diagonalization.","supporting_citations":[{"cited_title":"Computation o f unsteady nonlinear flows in cascades using a harmonic balance technique","cited_arxiv_id":null,"evidence_quote":"Supplies the harmonic balance form of the Navier-Stokes equations that the paper transforms and implements."},{"cited_title":"A new stabilized time-spectral finit e element solver for fast simulation of blood flow","cited_arxiv_id":null,"evidence_quote":"Prior frequency-domain stabilized finite element solver that this work extends, providing the GLS stabilization and pseudo-time-stepping scheme."},{"cited_title":"An augmented streamline upwind/p etrov-galerkin method for the time-spectral convection-diffusion equation","cited_arxiv_id":null,"evidence_quote":"Introduces the augmented stabilization method and analysis used to justify the GLS approach in the time-spectral setting."},{"cited_title":"The design and implementat ion of fftw3","cited_arxiv_id":null,"evidence_quote":"FFTW library implements the fast Fourier transforms that give the $O(N\\log N)$ matrix-vector products for the time-coupling matrix."},{"cited_title":"The vascular m odel repository: a pub- lic resource of medical imaging data and blood flow simulation results","cited_arxiv_id":null,"evidence_quote":"Supplies the patient-specific geometries used in the three test cases."},{"cited_title":"Troianowski, C","cited_arxiv_id":null,"evidence_quote":"Provides the Glenn procedure boundary conditions and clinically relevant simulation setup."},{"cited_title":"Patient-specific three-dimensional simulation of ldl accumulation in a human left coron ary artery in its healthy and atherosclerotic states","cited_arxiv_id":null,"evidence_quote":"Provides the coronary inflow profile and prior CFD setup for the left main coronary case."}],"review_version":1}