Journal of Computational Physics最新文献

{"title":"Algebraic dynamic multilevel (ADM) method for CO2 storage in heterogeneous saline aquifers","authors":"Mengjie Zhao ,&nbsp;Marc Gerritsma ,&nbsp;Mohammed Al Kobaisi ,&nbsp;Hadi Hajibeygi","doi":"10.1016/j.jcp.2025.114202","DOIUrl":"10.1016/j.jcp.2025.114202","url":null,"abstract":"<div><div>This work introduces a novel application of the Algebraic Dynamic Multilevel (ADM) method for simulating CO₂ storage in deep saline aquifers. By integrating a fully implicit coupling strategy, fully compositional thermodynamics, and adaptive mesh refinement, the ADM framework effectively models phenomena such as buoyancy-driven migration, convective dissolution, and phase partitioning under various subsurface conditions. The method starts with the construction of a hierarchy of multilevel grids and the generation of localized multiscale basis functions, which account for heterogeneities at each coarse level. During the simulation, the ADM method dynamically refines areas with significant overall CO₂ mass fraction gradients while coarsening smooth regions, thus optimizing computational resources without compromising the accuracy required to capture essential flow and transport characteristics. This dynamic grid adjustment is facilitated by algebraic prolongation and restriction operators, which map the fine-scale system onto a coarser grid suited to the evolving distribution of the CO₂ plume. This feature allows the ADM to navigate various coarsening thresholds efficiently, striking a trade-off between computational economy and detailed simulation accuracy. Even at relatively higher thresholds, key trapping mechanisms are captured with sufficient detail for quantification. These capabilities make the ADM framework well suited for long-term CO₂ sequestration in highly heterogeneous reservoirs, where large-scale models may otherwise become impractically expensive, offering a practical balance between the need for detailed simulations and manageable computational requirements. Overall, the ADM framework proves to be a robust tool for designing, monitoring, and analyzing large-scale CO₂ storage operations, supporting reliable and cost-effective solutions in carbon management.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114202"},"PeriodicalIF":3.8,"publicationDate":"2025-07-08","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144579476","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"Recursive sparse LU decomposition based on nested dissection and low rank approximations","authors":"Xuanru Zhu,&nbsp;Jun Lai","doi":"10.1016/j.jcp.2025.114231","DOIUrl":"10.1016/j.jcp.2025.114231","url":null,"abstract":"<div><div>When solving partial differential equations (PDEs) using finite difference or finite element methods, efficient solvers are required for handling large sparse linear systems. In this paper, a recursive sparse LU decomposition for matrices arising from the discretization of linear PDEs is proposed based on the nested dissection and low rank approximations. The matrix is reorganized based on the nested structure of the associated graph. After eliminating the interior vertices at the finest level, dense blocks on the separators are hierarchically sparsified using low rank approximations. To efficiently skeletonize these dense blocks, we split the separators into segments and introduce a hybrid algorithm to extract the low rank structures based on a randomized algorithm and the fast multipole method. The resulting decomposition yields a fast direct solver for sparse matrices, applicable to both symmetric and non-symmetric cases. Under a mild assumption on the compression rate of dense blocks, we prove an <span><math><mrow><mi>O</mi><mo>(</mo><mi>N</mi><mo>)</mo></mrow></math></span> complexity for the fast direct solver. Several numerical experiments are provided to verify the effectiveness of the proposed method.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114231"},"PeriodicalIF":3.8,"publicationDate":"2025-07-08","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144611753","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"A linearly implicit shock capturing scheme for compressible two-phase flows at all Mach numbers","authors":"Beatrice Battisti ,&nbsp;Walter Boscheri","doi":"10.1016/j.jcp.2025.114227","DOIUrl":"10.1016/j.jcp.2025.114227","url":null,"abstract":"<div><div>We present a semi-implicit solver for the solution of compressible two-phase flows governed by the Baer–Nunziato model. A novel linearly implicit discretization is proposed for the pressure fluxes as well as for the relaxation source terms, whereas an explicit scheme is retained for the nonlinear convective contributions. Consequently, the CFL-type stability condition on the maximum admissible time step is based only on the mean flow velocity and not on the sound speed of each phase, so that the novel scheme works uniformly for all Mach numbers. Central finite difference operators on Cartesian grids are adopted for the implicit terms, thus avoiding any need of numerical diffusion that might destroy accuracy in the low Mach number regime. To comply with high Mach number flows, shock capturing finite volume schemes are employed for the approximation of the convective fluxes. The discretization of the non-conservative terms ensures the preservation of moving equilibrium solutions, making the new method well-balanced. The new scheme is also proven to be asymptotic preserving in the low Mach limit of the mixture model. Second order of accuracy is achieved by means of an implicit-explicit (IMEX) time stepping algorithm combined with a total variation diminishing (TVD) reconstruction technique. The novel method is benchmarked against a set of test cases involving different Mach number regimes, permitting to validate both accuracy and robustness.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114227"},"PeriodicalIF":3.8,"publicationDate":"2025-07-08","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144605924","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"Solving continuum and rarefied flows using differentiable programming","authors":"Tianbai Xiao","doi":"10.1016/j.jcp.2025.114224","DOIUrl":"10.1016/j.jcp.2025.114224","url":null,"abstract":"<div><div>Accurate and efficient prediction of multi-scale flows remains a formidable challenge. Constructing theoretical models and numerical methods often involves the design and optimization of parameters. While gradient descent methods have been mainly manifested to shine in the wave of deep learning, composable automatic differentiation can advance scientific computing where the application of classical adjoint methods alone is infeasible or cumbersome. Differentiable programming provides a novel paradigm that unifies data structures and control flows and facilitates gradient-based optimization of parameters in a computer program. This paper addresses the notion and implementation of the first solution algorithm for multi-scale flow physics across continuum and rarefied regimes based on differentiable programming. The fully differentiable simulator provides a unified framework for the convergence of computational fluid dynamics and machine learning, i.e., scientific machine learning. Specifically, parameterized flow models and numerical methods can be constructed for forward physical processes, while the parameters can be trained on the fly with the help of the gradients that are taken through the backward passes of the whole simulation program, a.k.a., end-to-end optimization. As a result, versatile data-augmented modeling and simulation can be achieved for physics discovery, surrogate modeling, and simulation acceleration. The fundamentals and implementation of the solution algorithm are demonstrated in detail. Numerical experiments, including forward and inverse problems for hydrodynamic and kinetic equations, are presented to demonstrate the performance of the numerical method. The open-source codes to reproduce the numerical results are available under the MIT license.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114224"},"PeriodicalIF":3.8,"publicationDate":"2025-07-06","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144605923","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"High order treatment of moving curved boundaries: Arbitrary-Lagrangian-Eulerian methods with a shifted boundary polynomial correction","authors":"Walter Boscheri ,&nbsp;Mirco Ciallella","doi":"10.1016/j.jcp.2025.114215","DOIUrl":"10.1016/j.jcp.2025.114215","url":null,"abstract":"<div><div>In this paper we present a novel approach for the prescription of high order boundary conditions when approximating the solution of the Euler equations for compressible gas dynamics on curved moving domains. When dealing with curved boundaries, the consistency of boundary conditions is a real challenge, and it becomes even more challenging in the context of moving domains discretized with high order Arbitrary-Lagrangian-Eulerian (ALE) schemes. The ALE formulation is particularly well-suited for handling moving and deforming domains, thus allowing for the simulation of complex fluid-structure interaction problems. However, if not properly treated, the imposition of boundary conditions can lead to significant errors in the numerical solution, which can spoil the high order discretization of the underlying mathematical model. In order to tackle this issue, we propose a new method based on the recently developed shifted boundary polynomial correction, which was originally proposed in the discontinuous Galerkin (DG) framework on fixed meshes. The new method is integrated into the space-time corrector step of a direct ALE finite volume method to account for the local curvature of the moving boundary by only exploiting the high order reconstruction polynomial of the finite volume control volume. It relies on a correction based on the extrapolated value of the cell polynomial evaluated at the true geometry, thus not requiring the explicit evaluation of high order Taylor series. This greatly simplifies the treatment of moving curved boundaries, as it allows for the use of standard simplicial meshes, which are much easier to generate and move than curvilinear ones, especially for 3D time-dependent problems. Several numerical experiments are presented demonstrating the high order convergence properties of the new method in the context of compressible flows in moving curved domains, which remain approximated by piecewise linear elements.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114215"},"PeriodicalIF":3.8,"publicationDate":"2025-07-06","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144672526","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"Implementation of a Mesh refinement algorithm into the quasi-static PIC code QuickPIC","authors":"Q. Su ,&nbsp;F. Li ,&nbsp;W. An ,&nbsp;V. Decyk ,&nbsp;Y. Zhao ,&nbsp;L. Hildebrand ,&nbsp;T.N. Dalichaouch ,&nbsp;S. Zhou ,&nbsp;E.P. Alves ,&nbsp;A.S. Almgren ,&nbsp;W.B. Mori","doi":"10.1016/j.jcp.2025.114225","DOIUrl":"10.1016/j.jcp.2025.114225","url":null,"abstract":"<div><div>Plasma-based acceleration (PBA) has emerged as a promising candidate for the accelerator technology used to build a future linear collider and/or an advanced light source. In PBA, a trailing or witness particle beam is accelerated in the plasma wave wakefield (WF) created by a laser or particle beam driver. The WF is often nonlinear and involves the crossing of plasma particle trajectories in real space and thus particle-in-cell methods are used. The distance over which the drive beam evolves is several orders of magnitude larger than the wake wavelength. This large disparity in length scales is amenable to the quasi-static approach. Three-dimensional (3D), quasi-static (QS), particle-in-cell (PIC) codes, e.g., QuickPIC, have been shown to provide high fidelity simulation capability with 2-4 orders of magnitude speedup over 3D fully explicit PIC codes. In PBA, the witness beam needs to be matched to the focusing forces of the WF to reduce the emittance growth. In some linear collider designs, the matched spot size of the witness beam can be 2 to 3 orders of magnitude smaller than the spot size (and wavelength) of the wakefield. Such an additional disparity in length scales is ideal for mesh refinement where the WF within the witness beam is described on a finer mesh than the rest of the WF. A mesh refinement scheme is described that has been implemented into the 3D QS PIC code, QuickPIC. Very fine (high) resolution is used in a small spatial region that includes the witness beam and progressively coarser resolutions in the rest of the simulation domain. A fast multigrid Poisson solver has been implemented for the field solve on the refined meshes and a Fast Fourier Transform (FFT) based Poisson solver is used for the coarse mesh. The code has been parallelized with both MPI and OpenMP, and the parallel scalability has also been improved by using pipelining. A preliminary adaptive mesh refinement technique is described to optimize the computational time for simulations with an evolving witness beam size. Several test problems are used to verify that the mesh refinement algorithm provides accurate results. Additionally, the results are benchmarked against highly resolved simulations exhibiting near-azimuthal symmetry, performed using QPAD—a novel hybrid QS PIC code that uses a PIC description in the coordinates <span><math><mrow><mo>(</mo><mi>r</mi><mo>,</mo><mi>c</mi><mi>t</mi><mo>−</mo><mi>z</mi><mo>)</mo></mrow></math></span> and a gridless description in the azimuthal angle, <span><math><mi>ϕ</mi></math></span>.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114225"},"PeriodicalIF":3.8,"publicationDate":"2025-07-06","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144686017","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"The Helmholtz Dirichlet and Neumann problems on piecewise smooth open curves","authors":"Johan Helsing ,&nbsp;Shidong Jiang","doi":"10.1016/j.jcp.2025.114223","DOIUrl":"10.1016/j.jcp.2025.114223","url":null,"abstract":"<div><div>A numerical scheme is presented for solving the Helmholtz equation with Dirichlet or Neumann boundary conditions on piecewise smooth open curves, where the curves may have corners and multiple junctions. Existing integral equation methods for smooth open curves rely on analyzing the exact singularities of the density at endpoints for associated integral operators, explicitly extracting these singularities from the densities in the formulation, and using global quadrature to discretize the boundary integral equation. Extending these methods to handle curves with corners and multiple junctions is challenging because the singularity analysis becomes much more complex, and constructing high-order quadrature for discretizing layer potentials with singular and hypersingular kernels and singular densities is nontrivial. The proposed scheme is built upon the following two observations. First, the single-layer potential operator and the normal derivative of the double-layer potential operator serve as effective preconditioners for each other locally. Second, the recursively compressed inverse preconditioning (RCIP) method can be extended to address “implicit” second-kind integral equations. The scheme is high-order, adaptive, and capable of handling corners and multiple junctions without prior knowledge of the density singularity. It is also compatible with fast algorithms, such as the fast multipole method. The performance of the scheme is illustrated with several numerical examples.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114223"},"PeriodicalIF":3.8,"publicationDate":"2025-07-05","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144588491","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"High order discrete unified gas kinetic scheme with multi-moment constrained conservative semi-Lagrangian reconstruction","authors":"Jiaqi Bu ,&nbsp;Weidong Li ,&nbsp;Ming Fang ,&nbsp;Zhaoli Guo","doi":"10.1016/j.jcp.2025.114206","DOIUrl":"10.1016/j.jcp.2025.114206","url":null,"abstract":"<div><div>In this work, a high-order discrete unified gas kinetic (HDUGKS) scheme for multiscale gas flows was presented. The proposed scheme is a high-order finite volume method, which holds good conservation property of the scheme. Moreover, in the proposed scheme, the high-order accuracy was achieved by a constrained interpolation profile conservative semi-Lagrangian (CIP-CSL) reconstruction, where both point values (PVs) and volume-integrated averages (VIAs) are defined and used to fulfill a compact reconstruction of high order polynomials in a local finite volume cell. The PVs are updated using semi-Lagrangian scheme along the characteristics of the kinetic model equations to maintain high efficiency and low numerical dissipation, but the VIAs are evolved under the discrete unified gas kinetic scheme (DUGKS) framework to preserve the conservation. To demonstrate the idea of the HDUGKS, a cubic polynomial reconstruction based one-dimensional HDUGKS was presented. Furthermore, the high order accuracy property of the proposed scheme was evaluated by five one-dimensional numerical benchmarks including: (a) 1-D sine wave, (b) shock structure, (c) Sod’s shock tube problem, (d) Lax shock tube problem, (e) Shu-Osher problem, and the numerical results show that, while there is one order of accuracy reduction in the proposed scheme due to the low order semi-Lagrangian updating of the PVs, the proposed scheme can almost achieve third order accuracy and obtain more accurate results with fewer meshes than the DUGKS. Additionally, the numerical results indicate that the present method maintains the multiscale properties of the original DUGKS and serve as an efficient numerical method for flow simulations in all Knudsen number flow regime as well.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114206"},"PeriodicalIF":3.8,"publicationDate":"2025-07-04","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144572062","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"High order-accurate solution of scattering integral equations with unbounded solutions at corners","authors":"Constantine Sideris ,&nbsp;Davit Aslanyan ,&nbsp;Oscar P. Bruno","doi":"10.1016/j.jcp.2025.114213","DOIUrl":"10.1016/j.jcp.2025.114213","url":null,"abstract":"<div><div>Although high-order Maxwell integral equation solvers provide significant advantages in terms of speed and accuracy over corresponding low-order integral methods, their performance significantly degrades in presence of non-smooth geometries—owing to field enhancement and singularities that arise at sharp edges and corners which, if left untreated, give rise to significant accuracy losses. The problem is particularly challenging in cases where the density function—that is, the solution to the integral equation—exhibits unbounded singular behavior at edges and corners. While such difficulties can be circumvented in two-dimensional configurations, they constitute an intrinsic feature of existing three-dimensional Maxwell integral formulations, in which the tangential component of the surface current density diverges along edges. In order to tackle the problem this paper restricts attention to the simplest context in which the unbounded-density difficulty arises, namely, integral formulations in 2D space whose integral density blows up at corners; the strategies proposed, however, can be generalized to the 3D context. The novel methodologies presented in this paper yield high-order convergence for such challenging equations and achieve highly accurate solutions (even near edges and corners) without requiring a-priori analysis of the geometry or use of singular bases.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"539 ","pages":"Article 114213"},"PeriodicalIF":3.8,"publicationDate":"2025-07-04","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144605922","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}

{"title":"Multiscale approximation and two-grid preconditioner for extremely anisotropic heat flow","authors":"Maria Vasilyeva ,&nbsp;Golo Wimmer ,&nbsp;Ben S. Southworth","doi":"10.1016/j.jcp.2025.114201","DOIUrl":"10.1016/j.jcp.2025.114201","url":null,"abstract":"<div><div>We consider anis heat flow with extreme anisotropy, as arises in magnetized plasmas for fusion applications. Such problems pose significant challenges in both obtaining an accurate approximation as well in the construction of an efficient solver. In both cases, the underlying difficulty is in forming an accurate approximation of temperature fields that follow the direction of complex, non-grid-aligned magnetic fields. In this work, we construct a highly accurate coarse grid approximation using spectral multiscale basis functions based on local anisotropic normalized Laplacians. We show that the local generalized spectral problems yield local modes that align with magnetic fields, and provide an excellent coarse-grid approximation of the problem. We then utilize this spectral coarse space as an approximation in itself, and as the coarse-grid in a two-level spectral preconditioner. Numerical results are presented for several magnetic field distributions and anisotropy ratios up to <span><math><msup><mn>10</mn><mn>12</mn></msup></math></span>, showing highly accurate results with a large system size reduction, and two-grid preconditioning that converges in <span><math><mrow><mi>O</mi><mo>(</mo><mn>1</mn><mo>)</mo></mrow></math></span> iterations, independent of anisotropy. Aptara</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"538 ","pages":"Article 114201"},"PeriodicalIF":3.8,"publicationDate":"2025-07-02","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"144535934","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"物理与天体物理","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}