Pub Date : 2026-05-01Epub Date: 2025-12-23DOI: 10.1016/j.jcp.2025.114610
Tom Hickling , Jonathan F. MacArt , Justin Sirignano , Den Waidmann
<div><div>Turbulent flows are chaotic and unsteady, but their statistical distribution converges to a steady state. Engineering quantities of interest typically take the form of time-averaged statistics such as <span><math><mrow><mfrac><mn>1</mn><mi>t</mi></mfrac><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup><mi>f</mi><mrow><mo>(</mo><mi>u</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>τ</mi><mo>;</mo><mi>θ</mi><mo>)</mo></mrow><mo>)</mo></mrow><mi>d</mi><mi>τ</mi><mover><mo>→</mo><mrow><mi>t</mi><mo>→</mo><mi>∞</mi></mrow></mover><mi>F</mi><mrow><mo>(</mo><mi>x</mi><mo>;</mo><mi>θ</mi><mo>)</mo></mrow></mrow></math></span>, where <em>u</em>(<em>x, t</em>; <em>θ</em>) is a solution of the Navier-Stokes equations with parameters <em>θ</em>. Optimizing over the time-averaged statistic <em>F</em>(<em>x</em>; <em>θ</em>) has many engineering applications, including geometric optimization, flow control, and closure modeling. However, optimizing <em>F</em>(<em>x</em>; <em>θ</em>) is non-trivial and remains an open challenge. The fundamental obstacle is the chaoticity of turbulent flows: gradients calculated with the adjoint method diverge exponentially as <em>t</em> → ∞, and existing stabilized approaches (such as least squares shadowing) are incapable of scaling to physically representative mesh resolutions, which can require over 10<sup>7</sup> degrees of freedom (the number of PDE variables × mesh points for finite-difference methods).</div><div>We develop a new online gradient-flow (OGF) method that enables optimizing over the steady-state statistics of chaotic, unsteady, turbulence-resolving simulations. The method can be viewed as a stochastic gradient descent method for chaotic systems, where a noisy <em>online</em> estimate is forward-propagated for the gradient of <em>F</em>(<em>x</em>; <em>θ</em>), and online updates of the parameters <em>θ</em> are performed simultaneously. The fully online nature of the algorithm facilitates faster optimization progress, and its combination with a finite-difference estimate avoids diverging gradients. The cost of the method is linear in both the degrees of freedom and the number of parameters. Although the cost of existing stabilized adjoint approaches is independent of the number of parameters, they are incapable of scaling to high degree of freedom systems. OGF therefore fills a critical capability gap: optimization of intermediate numbers of parameters over large degree of freedom chaotic systems.</div><div>We demonstrate the OGF method in optimizations over four different chaotic systems: the Lorenz-63 equation, including training an ODE-embedded neural network with 501 parameters; the Kuramoto-Sivashinsky equation; Navier-Stokes solutions of compressible, forced homogeneous isotropic turbulence; and a seven-parameter perceptron-based subgrid-scale model in large eddy simulations of turbulent channel flow. The latter two systems have 10<sup>7</sup>–10<sup>8</sup> degrees of freedom, showing the scalab
{"title":"OGF: An online gradient flow method for optimizing the statistical steady-state time averages of unsteady turbulent flows","authors":"Tom Hickling , Jonathan F. MacArt , Justin Sirignano , Den Waidmann","doi":"10.1016/j.jcp.2025.114610","DOIUrl":"10.1016/j.jcp.2025.114610","url":null,"abstract":"<div><div>Turbulent flows are chaotic and unsteady, but their statistical distribution converges to a steady state. Engineering quantities of interest typically take the form of time-averaged statistics such as <span><math><mrow><mfrac><mn>1</mn><mi>t</mi></mfrac><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup><mi>f</mi><mrow><mo>(</mo><mi>u</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>τ</mi><mo>;</mo><mi>θ</mi><mo>)</mo></mrow><mo>)</mo></mrow><mi>d</mi><mi>τ</mi><mover><mo>→</mo><mrow><mi>t</mi><mo>→</mo><mi>∞</mi></mrow></mover><mi>F</mi><mrow><mo>(</mo><mi>x</mi><mo>;</mo><mi>θ</mi><mo>)</mo></mrow></mrow></math></span>, where <em>u</em>(<em>x, t</em>; <em>θ</em>) is a solution of the Navier-Stokes equations with parameters <em>θ</em>. Optimizing over the time-averaged statistic <em>F</em>(<em>x</em>; <em>θ</em>) has many engineering applications, including geometric optimization, flow control, and closure modeling. However, optimizing <em>F</em>(<em>x</em>; <em>θ</em>) is non-trivial and remains an open challenge. The fundamental obstacle is the chaoticity of turbulent flows: gradients calculated with the adjoint method diverge exponentially as <em>t</em> → ∞, and existing stabilized approaches (such as least squares shadowing) are incapable of scaling to physically representative mesh resolutions, which can require over 10<sup>7</sup> degrees of freedom (the number of PDE variables × mesh points for finite-difference methods).</div><div>We develop a new online gradient-flow (OGF) method that enables optimizing over the steady-state statistics of chaotic, unsteady, turbulence-resolving simulations. The method can be viewed as a stochastic gradient descent method for chaotic systems, where a noisy <em>online</em> estimate is forward-propagated for the gradient of <em>F</em>(<em>x</em>; <em>θ</em>), and online updates of the parameters <em>θ</em> are performed simultaneously. The fully online nature of the algorithm facilitates faster optimization progress, and its combination with a finite-difference estimate avoids diverging gradients. The cost of the method is linear in both the degrees of freedom and the number of parameters. Although the cost of existing stabilized adjoint approaches is independent of the number of parameters, they are incapable of scaling to high degree of freedom systems. OGF therefore fills a critical capability gap: optimization of intermediate numbers of parameters over large degree of freedom chaotic systems.</div><div>We demonstrate the OGF method in optimizations over four different chaotic systems: the Lorenz-63 equation, including training an ODE-embedded neural network with 501 parameters; the Kuramoto-Sivashinsky equation; Navier-Stokes solutions of compressible, forced homogeneous isotropic turbulence; and a seven-parameter perceptron-based subgrid-scale model in large eddy simulations of turbulent channel flow. The latter two systems have 10<sup>7</sup>–10<sup>8</sup> degrees of freedom, showing the scalab","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114610"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146170754","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}
Pub Date : 2026-05-01Epub Date: 2026-01-22DOI: 10.1016/j.jcp.2026.114701
A. Colaïtis , S. Guisset , J. Breil
Many applications of physics and engineering involve wide ranges of time and spatial scales. The numerical simulation of localized small scales such as shock waves and material interfaces requires a large number of computational cells in these regions. For these applications, Lagrangian and Arbitrary-Lagrangian-Eulerian (ALE) related methods are engaging since the moving mesh feature naturally brings mesh cells on shock discontinuities and material interfaces are carefully captured. In addition, Adaptive-Mesh-Refinement (AMR) strategies aim to optimize computational resources by concentrating finer mesh cells only in areas of interest while using coarser cells elsewhere. A key but challenging AMR requirement consists in efficiently distributing the computational effort to achieve high accuracy without the prohibitive computational costs associated with uniformly fine grids. In this document, the coupling of the p4est AMR library with a cell-centered Lagrangian scheme is presented with the goal to perform reliable 3D Lagrangian-AMR and indirect Euler-AMR multi-material simulations. In particular, it is shown that starting from a 3D indirect ALE code, the memory management and load balancing requirements can be delegated to an external library (here the p4est library) to unlock ALE-AMR capabilities. First, we present a strategy to transcribe the octant-based connectivity of the 3D AMR framework with that of an unstructured mesh of polygonal cells used in Lagrangian hydrodynamics. Then, we show how refinement and coarsening operations must be adapted to the particular Lagrangian framework to ensure the conservation of volume during those steps. Finally, several numerical test cases are presented that demonstrate the capabilities of the Lagrangian-AMR and indirect Euler-AMR algorithms.
{"title":"A cell-centered AMR-ALE framework for 3D multi-material hydrodynamics. Part I: Lagrangian and indirect Euler AMR algorithms","authors":"A. Colaïtis , S. Guisset , J. Breil","doi":"10.1016/j.jcp.2026.114701","DOIUrl":"10.1016/j.jcp.2026.114701","url":null,"abstract":"<div><div>Many applications of physics and engineering involve wide ranges of time and spatial scales. The numerical simulation of localized small scales such as shock waves and material interfaces requires a large number of computational cells in these regions. For these applications, Lagrangian and Arbitrary-Lagrangian-Eulerian (ALE) related methods are engaging since the moving mesh feature naturally brings mesh cells on shock discontinuities and material interfaces are carefully captured. In addition, Adaptive-Mesh-Refinement (AMR) strategies aim to optimize computational resources by concentrating finer mesh cells only in areas of interest while using coarser cells elsewhere. A key but challenging AMR requirement consists in efficiently distributing the computational effort to achieve high accuracy without the prohibitive computational costs associated with uniformly fine grids. In this document, the coupling of the <em>p4est</em> AMR library with a cell-centered Lagrangian scheme is presented with the goal to perform reliable 3D Lagrangian-AMR and indirect Euler-AMR multi-material simulations. In particular, it is shown that starting from a 3D indirect ALE code, the memory management and load balancing requirements can be delegated to an external library (here the <em>p4est</em> library) to unlock ALE-AMR capabilities. First, we present a strategy to transcribe the octant-based connectivity of the 3D AMR framework with that of an unstructured mesh of polygonal cells used in Lagrangian hydrodynamics. Then, we show how refinement and coarsening operations must be adapted to the particular Lagrangian framework to ensure the conservation of volume during those steps. Finally, several numerical test cases are presented that demonstrate the capabilities of the Lagrangian-AMR and indirect Euler-AMR algorithms.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114701"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146075600","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}
Pub Date : 2026-05-01Epub Date: 2026-01-24DOI: 10.1016/j.jcp.2026.114705
Jingfei Chen , Minxin Chen , Jingrun Chen
Transformer performs remarkably across a diverse array of natural language processing tasks. Its applicability has been extended to the domain of partial differential equations, giving rise to two novel models: Fourier and Galerkin. However, the self-attention module inherent in Transformers exhibits a quadratic computational complexity with respect to the input sequence length n, as observed in the Fourier model, leading to substantial computational overheads for long input sequences. A detailed analysis illustrates that the correlation matrix in the self-attention mechanism exhibits a low-rank property. We incorporate this structure into the network architecture, where the self-attention mechanism degenerates into an adaptive low-rank ResNet-type network (ALRN). This network can adaptively capture the rank of the correlation matrix. As a consequence, the ALRN model entails a computational complexity of in contrast to for the Fourier model and for the Galerkin model, where k denotes the rank of the correlation matrix and d represents the dimension of the feature space. Concerning parameter space, it is noteworthy that Fourier and Galerkin methods necessitate , while the ALRN model demands . Therefore, the ALRN model offers a clear advantage when . Numerical results for Burgers’ equation, Darcy flow, the inverse coefficient identification for the Darcy flow, and the Navier-Stokes equation demonstrate superior efficiency and require fewer parameters while maintaining accuracy.
{"title":"An operator learning method for solving partial differential equations: From transformer to adaptive low-rank resnet-type network","authors":"Jingfei Chen , Minxin Chen , Jingrun Chen","doi":"10.1016/j.jcp.2026.114705","DOIUrl":"10.1016/j.jcp.2026.114705","url":null,"abstract":"<div><div>Transformer performs remarkably across a diverse array of natural language processing tasks. Its applicability has been extended to the domain of partial differential equations, giving rise to two novel models: Fourier and Galerkin. However, the self-attention module inherent in Transformers exhibits a quadratic computational complexity with respect to the input sequence length <em>n</em>, as observed in the Fourier model, leading to substantial computational overheads for long input sequences. A detailed analysis illustrates that the correlation matrix in the self-attention mechanism exhibits a low-rank property. We incorporate this structure into the network architecture, where the self-attention mechanism degenerates into an adaptive low-rank ResNet-type network (ALRN). This network can adaptively capture the rank of the correlation matrix. As a consequence, the ALRN model entails a computational complexity of <span><math><mrow><mi>O</mi><mo>(</mo><mi>n</mi><mrow><mo>(</mo><mn>2</mn><msup><mi>d</mi><mn>2</mn></msup><mo>+</mo><mn>4</mn><mi>k</mi><mi>d</mi><mo>)</mo></mrow><mo>)</mo></mrow></math></span> in contrast to <span><math><mrow><mi>O</mi><mo>(</mo><mn>4</mn><msup><mi>n</mi><mn>2</mn></msup><mi>d</mi><mo>)</mo></mrow></math></span> for the Fourier model and <span><math><mrow><mi>O</mi><mo>(</mo><mn>12</mn><mi>n</mi><msup><mi>d</mi><mn>2</mn></msup><mo>)</mo></mrow></math></span> for the Galerkin model, where <em>k</em> denotes the rank of the correlation matrix and <em>d</em> represents the dimension of the feature space. Concerning parameter space, it is noteworthy that Fourier and Galerkin methods necessitate <span><math><mrow><mi>O</mi><mo>(</mo><mn>3</mn><msup><mi>d</mi><mn>2</mn></msup><mo>)</mo></mrow></math></span>, while the ALRN model demands <span><math><mrow><mi>O</mi><mo>(</mo><msup><mi>d</mi><mn>2</mn></msup><mo>+</mo><mn>2</mn><mi>n</mi><mi>k</mi><mo>)</mo></mrow></math></span>. Therefore, the ALRN model offers a clear advantage when <span><math><mrow><mi>n</mi><mo><</mo><mi>O</mi><mo>(</mo><msup><mi>d</mi><mn>2</mn></msup><mo>)</mo></mrow></math></span>. Numerical results for Burgers’ equation, Darcy flow, the inverse coefficient identification for the Darcy flow, and the Navier-Stokes equation demonstrate superior efficiency and require fewer parameters while maintaining accuracy.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114705"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146075599","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}
Pub Date : 2026-05-01Epub Date: 2026-01-22DOI: 10.1016/j.jcp.2026.114697
Spencer Lee , Daniel Appelo
This work introduces the High-Order Hermite Optimization (HOHO) method, an open-loop discrete adjoint method for quantum optimal control. Our method is the first of its kind to efficiently compute exact (discrete) gradients when using continuous, parameterized control pulses while solving the forward equations (e.g. Schrodinger’s equation or the Linblad master equation) with an arbitrarily high-order Hermite Runge-Kutta method. The HOHO method is implemented in QuantumGateDesign.jl, an open-source software package for the Julia programming language, which we use to perform numerical experiments comparing the method to Juqbox.jl. For realistic model problems we observe speedups up to 775x.
{"title":"High-order Hermite optimization: Fast and exact gradient computation in open-loop quantum optimal control using a discrete adjoint approach","authors":"Spencer Lee , Daniel Appelo","doi":"10.1016/j.jcp.2026.114697","DOIUrl":"10.1016/j.jcp.2026.114697","url":null,"abstract":"<div><div>This work introduces the <em>High-Order Hermite Optimization</em> (HOHO) method, an open-loop discrete adjoint method for quantum optimal control. Our method is the first of its kind to efficiently compute exact (discrete) gradients when using continuous, parameterized control pulses while solving the forward equations (e.g. Schrodinger’s equation or the Linblad master equation) with an <em>arbitrarily high-order</em> Hermite Runge-Kutta method. The HOHO method is implemented in <em>QuantumGateDesign.jl</em>, an open-source software package for the Julia programming language, which we use to perform numerical experiments comparing the method to Juqbox.jl. For realistic model problems we observe speedups up to 775x.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114697"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146075598","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}
Pub Date : 2026-05-01Epub Date: 2026-01-22DOI: 10.1016/j.jcp.2026.114699
Chang Wei , Haotian Liu , Shangming Li
The immersed finite element/difference (IFED) method is a numerical framework for modeling fluid-structure interaction in single-phase flows. In this work, the IFED method is extended to multiphase FSI with emphasis on deformable structures. The proposed approach preserves the geometric flexibility for deformable structures and integrates the level set method to capture phase interfaces. The gas-liquid interface is tracked by advecting a level set, whereas the structural level set is explicitly reconstructed to produce a discrete signed distance field on the Cartesian grid. Once the level set fields are available, fluid density and viscosity are assigned via regularized Heaviside functions to smooth interfacial transitions. Fluid-structure coupling follows the IFED method, in which a force-spreading operator spreads Lagrangian forces onto the Cartesian grid, and a velocity-restriction operator transfers Eulerian velocities back to the structural mesh. These operators are implemented within the multiphase framework. A limitation of the original IFED formulation is the homogeneous time-step coupling that imposes the same time step on fluid and solid subdomains. To address this, a time-splitting scheme is proposed using two Lagrangian representations: the original mesh interacts with the fluid without structural constitutive response, and the auxiliary mesh is coupled to the former via a penalty force formulation. This allows the structural elastodynamic equations to be solved multiple times within each fluid time step using a standard Galerkin finite element method on the auxiliary mesh. The resulting material response is transmitted back to the original mesh as another penalty body force, which is subsequently used to compute the Eulerian force density. The proposed scheme is applicable to both the original IFED method and its multiphase extension. Two 2D dam-break tests indicate that the multiphase IFED method provides accurate predictions for deformable structures with low to moderate stiffness. For high-stiffness structures, the time-splitting scheme achieves achieves a speedup of about 3.8 times in the Turek-Hron test at Pa relative to the original version, solely by increasing structural substeps. A dam-break impact test involving a deformable body with Pa further demonstrates the effectiveness of the time-splitting scheme. The study provides new insights into the modeling of multiphase FSI problems involving deformable structures.
{"title":"A multiphase IFED method for fluid-structure interactions","authors":"Chang Wei , Haotian Liu , Shangming Li","doi":"10.1016/j.jcp.2026.114699","DOIUrl":"10.1016/j.jcp.2026.114699","url":null,"abstract":"<div><div>The immersed finite element/difference (IFED) method is a numerical framework for modeling fluid-structure interaction in single-phase flows. In this work, the IFED method is extended to multiphase FSI with emphasis on deformable structures. The proposed approach preserves the geometric flexibility for deformable structures and integrates the level set method to capture phase interfaces. The gas-liquid interface is tracked by advecting a level set, whereas the structural level set is explicitly reconstructed to produce a discrete signed distance field on the Cartesian grid. Once the level set fields are available, fluid density and viscosity are assigned via regularized Heaviside functions to smooth interfacial transitions. Fluid-structure coupling follows the IFED method, in which a force-spreading operator spreads Lagrangian forces onto the Cartesian grid, and a velocity-restriction operator transfers Eulerian velocities back to the structural mesh. These operators are implemented within the multiphase framework. A limitation of the original IFED formulation is the homogeneous time-step coupling that imposes the same time step on fluid and solid subdomains. To address this, a time-splitting scheme is proposed using two Lagrangian representations: the original mesh interacts with the fluid without structural constitutive response, and the auxiliary mesh is coupled to the former via a penalty force formulation. This allows the structural elastodynamic equations to be solved multiple times within each fluid time step using a standard Galerkin finite element method on the auxiliary mesh. The resulting material response is transmitted back to the original mesh as another penalty body force, which is subsequently used to compute the Eulerian force density. The proposed scheme is applicable to both the original IFED method and its multiphase extension. Two 2D dam-break tests indicate that the multiphase IFED method provides accurate predictions for deformable structures with low to moderate stiffness. For high-stiffness structures, the time-splitting scheme achieves achieves a speedup of about 3.8 times in the Turek-Hron test at <span><math><mrow><msub><mi>G</mi><mi>s</mi></msub><mo>=</mo><mn>2</mn><mo>×</mo><msup><mn>10</mn><mn>7</mn></msup></mrow></math></span> Pa relative to the original version, solely by increasing structural substeps. A dam-break impact test involving a deformable body with <span><math><mrow><msub><mi>E</mi><mi>s</mi></msub><mo>=</mo><mn>5.0</mn><mo>×</mo><msup><mn>10</mn><mn>10</mn></msup></mrow></math></span> Pa further demonstrates the effectiveness of the time-splitting scheme. The study provides new insights into the modeling of multiphase FSI problems involving deformable structures.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114699"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146075597","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}
The simulation of flows presenting contact discontinuities, vorticity, and large variations in spatial scales can be performed in a framework coupling Arbitrary Lagrangian Eulerian (ALE) algorithms and Adaptive Mesh Refinement (AMR). This coupling requires adaptation of ALE rezoning techniques to meshes containing nonconformal nodes arising from both the AMR topology and the junction of mesh blocks. In this paper, we present an ALE rezoning strategy that is compatible with such meshes, and that can also act as a disentangling algorithm. Emphasis is put on an algorithm that respects intrinsic Lagrangian mesh properties in order to preserve accuracy around discontinuities. To that end, we adapt the weighted linesweep algorithm to nonconformal block-structured AMR meshes. Then, we present control parameters introduced in the method for it to be applicable in practical situations. Notably, the method is coupled to a specific metric optimization in order to palliate some shortcomings of the linesweep method. Finally, numerical test cases are presented that feature the capabilities of the ALE-AMR algorithm for flows that present discontinuities, vorticity, and a variety of scales. Notably, we show that our ALE-AMR algorithm gives results at least similar to Euler-AMR, but provides better accuracy in cases where discontinuities are involved, thanks to a method that respects the Lagrangian features of the mesh. Additionally, it enables Euler-AMR-like computations on domains with temporally varying domain boundaries.
{"title":"A cell-centered AMR-ALE framework for 3D multi-material hydrodynamics. Part II: linesweep ALE rezoning for nonconformal block-structured AMR meshes","authors":"Arnaud Colaïtis , Sébastien Guisset , Jérôme Breil","doi":"10.1016/j.jcp.2026.114702","DOIUrl":"10.1016/j.jcp.2026.114702","url":null,"abstract":"<div><div>The simulation of flows presenting contact discontinuities, vorticity, and large variations in spatial scales can be performed in a framework coupling Arbitrary Lagrangian Eulerian (ALE) algorithms and Adaptive Mesh Refinement (AMR). This coupling requires adaptation of ALE rezoning techniques to meshes containing nonconformal nodes arising from both the AMR topology and the junction of mesh blocks. In this paper, we present an ALE rezoning strategy that is compatible with such meshes, and that can also act as a disentangling algorithm. Emphasis is put on an algorithm that respects intrinsic Lagrangian mesh properties in order to preserve accuracy around discontinuities. To that end, we adapt the weighted linesweep algorithm to nonconformal block-structured AMR meshes. Then, we present control parameters introduced in the method for it to be applicable in practical situations. Notably, the method is coupled to a specific metric optimization in order to palliate some shortcomings of the linesweep method. Finally, numerical test cases are presented that feature the capabilities of the ALE-AMR algorithm for flows that present discontinuities, vorticity, and a variety of scales. Notably, we show that our ALE-AMR algorithm gives results at least similar to Euler-AMR, but provides better accuracy in cases where discontinuities are involved, thanks to a method that respects the Lagrangian features of the mesh. Additionally, it enables Euler-AMR-like computations on domains with temporally varying domain boundaries.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114702"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146025850","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}
Pub Date : 2026-05-01Epub Date: 2026-01-25DOI: 10.1016/j.jcp.2026.114706
Ignacio Puiggros T․ , A. Srikantha Phani
There has been increasing interest in methodologies that incorporate physics priors into neural network architectures to enhance their modeling capabilities. A family of these methodologies that has gained traction are Hamiltonian neural networks (HNN) and their variants. These architectures explicitly encode Hamiltonian mechanics both in their structure and loss function. Although Hamiltonian systems under nonholonomic constraints are in general not Hamiltonian, it is possible to formulate them in pseudo-Hamiltonian form, equipped with a Lie bracket which is almost Poisson. This opens the possibility of using some principles of HNNs in systems under nonholonomic constraints. The goal of the present work is to develop a modified Hamiltonian neural network architecture capable of modeling Hamiltonian systems under holonomic and nonholonomic constraints. A three-network parallel architecture is proposed to simultaneously learn the Hamiltonian of the system, the constraints, and their associated multipliers. A rolling disk and a ball on a spinning table are considered as canonical examples to assess the performance of the proposed Hamiltonian architecture. The numerical experiments are then repeated with a noisy training set to study modeling performance under more realistic conditions.
{"title":"Hamiltonian-based neural networks for systems under nonholonomic constraints","authors":"Ignacio Puiggros T․ , A. Srikantha Phani","doi":"10.1016/j.jcp.2026.114706","DOIUrl":"10.1016/j.jcp.2026.114706","url":null,"abstract":"<div><div>There has been increasing interest in methodologies that incorporate physics priors into neural network architectures to enhance their modeling capabilities. A family of these methodologies that has gained traction are Hamiltonian neural networks (HNN) and their variants. These architectures explicitly encode Hamiltonian mechanics both in their structure and loss function. Although Hamiltonian systems under nonholonomic constraints are in general not Hamiltonian, it is possible to formulate them in pseudo-Hamiltonian form, equipped with a Lie bracket which is almost Poisson. This opens the possibility of using some principles of HNNs in systems under nonholonomic constraints. The goal of the present work is to develop a modified Hamiltonian neural network architecture capable of modeling Hamiltonian systems under holonomic and nonholonomic constraints. A three-network parallel architecture is proposed to simultaneously learn the Hamiltonian of the system, the constraints, and their associated multipliers. A rolling disk and a ball on a spinning table are considered as canonical examples to assess the performance of the proposed Hamiltonian architecture. The numerical experiments are then repeated with a noisy training set to study modeling performance under more realistic conditions.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114706"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146170750","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}
Pub Date : 2026-05-01Epub Date: 2026-01-26DOI: 10.1016/j.jcp.2026.114708
Peiran Zhang
We consider the numerical reconstruction of the spatially dependent conductivity coefficient and the source term in elliptic partial differential equations in a two-dimensional convex polygonal domain, with homogeneous Dirichlet boundary condition and given interior observations of the solution. Using continuous data assimilation, we derive approximated gradients of the error function to update the reconstructed coefficients, which, in particular, avoids solving adjoint problems. New L2 error estimates are provided for the spatially discretized reconstructions. Numerical examples are given to illustrate the effectiveness of the method and demonstrate the error estimates. The numerical results also reveal a notable feature that the reconstruction is very robust to errors in coefficients.
{"title":"Numerical reconstruction of coefficients in elliptic equations using continuous data assimilation","authors":"Peiran Zhang","doi":"10.1016/j.jcp.2026.114708","DOIUrl":"10.1016/j.jcp.2026.114708","url":null,"abstract":"<div><div>We consider the numerical reconstruction of the spatially dependent conductivity coefficient and the source term in elliptic partial differential equations in a two-dimensional convex polygonal domain, with homogeneous Dirichlet boundary condition and given interior observations of the solution. Using continuous data assimilation, we derive approximated gradients of the error function to update the reconstructed coefficients, which, in particular, avoids solving adjoint problems. New <em>L</em><sup>2</sup> error estimates are provided for the spatially discretized reconstructions. Numerical examples are given to illustrate the effectiveness of the method and demonstrate the error estimates. The numerical results also reveal a notable feature that the reconstruction is very robust to errors in coefficients.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114708"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146075596","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}
Pub Date : 2026-05-01Epub Date: 2026-01-14DOI: 10.1016/j.jcp.2026.114683
James F. Kelly , Felipe A. V. De Bragança Alves , Stephen D. Eckermann , Francis X. Giraldo , P. Alex Reinecke , John T. Emmert
This paper presents and tests a deep-atmosphere, nonhydrostatic dynamical core (DyCore) targeted towards ground-to-thermosphere atmospheric prediction using the spectral element method (SEM) with IMplicit-EXplicit (IMEX) and Horizontally Explicit Vertically Implicit (HEVI) time-integration. Two versions of the DyCore are discussed, each based on a different formulation of the specific internal energy and continuity equations, which, unlike the dynamical cores developed for low-altitude atmospheric applications, are valid for variable composition atmospheres. The first version, which uses a product-rule (PR) form of the continuity and specific internal energy equations, contains an additional pressure dilation term and does not conserve mass. The second version, which does not use the product-rule (no-PR) in the continuity and specific internal energy, contains two terms to represent pressure dilation in the energy equation and conserves mass to machine precision regardless of time truncation error. The pressure gradient and gravitational forces in the momentum balance equation are reformulated to reduce numerical errors at high altitudes. These new equation sets were implemented in two SEM-based atmospheric models: the Nonhydrostatic Unified Model of the Atmosphere (NUMA) and the Navy Environmental Prediction sysTem Utilizing a Nonhydrostatic Engine (NEPTUNE). Numerical results using both a deep-atmosphere and shallow-atmosphere baroclinic instability, a balanced zonal flow, and a high-altitude orographic gravity wave verify the fidelity of the dynamics at low and high altitudes and for constant and variable composition atmospheres. These results are compared to those from existing deep-atmosphere dynamical cores and a Fourier-ray code, indicating that the proposed discretized equation sets are viable DyCore candidates for next-generation ground-to-thermosphere atmospheric models.
{"title":"A nonhydrostatic mass-conserving dynamical core for deep atmospheres of variable composition","authors":"James F. Kelly , Felipe A. V. De Bragança Alves , Stephen D. Eckermann , Francis X. Giraldo , P. Alex Reinecke , John T. Emmert","doi":"10.1016/j.jcp.2026.114683","DOIUrl":"10.1016/j.jcp.2026.114683","url":null,"abstract":"<div><div>This paper presents and tests a deep-atmosphere, nonhydrostatic dynamical core (DyCore) targeted towards ground-to-thermosphere atmospheric prediction using the spectral element method (SEM) with IMplicit-EXplicit (IMEX) and Horizontally Explicit Vertically Implicit (HEVI) time-integration. Two versions of the DyCore are discussed, each based on a different formulation of the specific internal energy and continuity equations, which, unlike the dynamical cores developed for low-altitude atmospheric applications, are valid for variable composition atmospheres. The first version, which uses a product-rule (PR) form of the continuity and specific internal energy equations, contains an additional pressure dilation term and does not conserve mass. The second version, which does not use the product-rule (no-PR) in the continuity and specific internal energy, contains two terms to represent pressure dilation in the energy equation and conserves mass to machine precision regardless of time truncation error. The pressure gradient and gravitational forces in the momentum balance equation are reformulated to reduce numerical errors at high altitudes. These new equation sets were implemented in two SEM-based atmospheric models: the Nonhydrostatic Unified Model of the Atmosphere (NUMA) and the Navy Environmental Prediction sysTem Utilizing a Nonhydrostatic Engine (NEPTUNE). Numerical results using both a deep-atmosphere and shallow-atmosphere baroclinic instability, a balanced zonal flow, and a high-altitude orographic gravity wave verify the fidelity of the dynamics at low and high altitudes and for constant and variable composition atmospheres. These results are compared to those from existing deep-atmosphere dynamical cores and a Fourier-ray code, indicating that the proposed discretized equation sets are viable DyCore candidates for next-generation ground-to-thermosphere atmospheric models.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114683"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146025858","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}
Pub Date : 2026-05-01Epub Date: 2026-01-20DOI: 10.1016/j.jcp.2026.114695
Cagatay Guventurk, Mehmet Sahin
An arbitrary Lagrangian Eulerian (ALE) framework presented in A mass conserving arbitrary Lagrangian-Eulerian formulation for three-dimensional multiphase fluid flows, International Journal for Numerical Methods in Fluids 94 (4), 346–376 has been extended to solve incompressible multiphase viscoelastic flow problems in two- and three-dimensions. The incompressible, isothermal linear momentum balance equations, coupled with the viscoelastic constitutive models Oldroyd-B and FENE-CR, are discretized using a div-stable, side-centered finite volume approach in which the velocity components are defined at the mid-points of element faces, the displacement vector is defined at the vertices, and the pressure and modified conformation tensor are defined at the element centroids. At the interface, the surface tension force is treated as a force tangential to the interface, and its normal vector is evaluated by using the mean weighted by sine and edge length reciprocals (MWSELR) approach. In order to ensure mass conservation of both species at machine precision, special attention is given to enforcing the kinematic boundary condition at the interface in the normal direction, while obeying the discrete geometric conservation law (DGCL). The numerical approach allows discontinuities in material properties, including density and viscosity, as well as in the pressure and modified conformation tensor across the interface. The discrete algebraic equations arising from the incompressible linear momentum balance equations are solved monolithically using a block preconditioner based on the BoomerAMG parallel algebraic multigrid solver from the HYPRE library, interfaced through PETSc. To validate the numerical algorithm, the benchmark problem of a single Newtonian or viscoelastic bubble (modeled using Oldroyd-B and FENE-CR) rising through a quiescent Newtonian or viscoelastic fluid is examined in both two- and three-dimensions. The numerical simulations exhibit excellent agreement with previous results in the literature and show strong consistency with mesh refinement. Positive and negative transient wakes are observed behind the bubble, demonstrating that the formation of a transient negative wake does not require a viscoelastic fluid model with shear-thinning behavior. The numerical approach successfully preserves the volume of the bubble to nearly machine precision and accurately captures discontinuities in the pressure and modified conformation tensor across the interface, where there are jumps in density and viscosity.
{"title":"An exact mass-conserving arbitrary Lagrangian-Eulerian framework for viscoelastic multiphase fluid flows","authors":"Cagatay Guventurk, Mehmet Sahin","doi":"10.1016/j.jcp.2026.114695","DOIUrl":"10.1016/j.jcp.2026.114695","url":null,"abstract":"<div><div>An arbitrary Lagrangian Eulerian (ALE) framework presented in <em>A mass conserving arbitrary Lagrangian-Eulerian formulation for three-dimensional multiphase fluid flows</em>, International Journal for Numerical Methods in Fluids 94 (4), 346–376 has been extended to solve incompressible multiphase viscoelastic flow problems in two- and three-dimensions. The incompressible, isothermal linear momentum balance equations, coupled with the viscoelastic constitutive models Oldroyd-B and FENE-CR, are discretized using a div-stable, side-centered finite volume approach in which the velocity components are defined at the mid-points of element faces, the displacement vector is defined at the vertices, and the pressure and modified conformation tensor are defined at the element centroids. At the interface, the surface tension force is treated as a force tangential to the interface, and its normal vector is evaluated by using the mean weighted by sine and edge length reciprocals (MWSELR) approach. In order to ensure mass conservation of both species at machine precision, special attention is given to enforcing the kinematic boundary condition at the interface in the normal direction, while obeying the discrete geometric conservation law (DGCL). The numerical approach allows discontinuities in material properties, including density and viscosity, as well as in the pressure and modified conformation tensor across the interface. The discrete algebraic equations arising from the incompressible linear momentum balance equations are solved monolithically using a block preconditioner based on the BoomerAMG parallel algebraic multigrid solver from the HYPRE library, interfaced through PETSc. To validate the numerical algorithm, the benchmark problem of a single Newtonian or viscoelastic bubble (modeled using Oldroyd-B and FENE-CR) rising through a quiescent Newtonian or viscoelastic fluid is examined in both two- and three-dimensions. The numerical simulations exhibit excellent agreement with previous results in the literature and show strong consistency with mesh refinement. Positive and negative transient wakes are observed behind the bubble, demonstrating that the formation of a transient negative wake does not require a viscoelastic fluid model with shear-thinning behavior. The numerical approach successfully preserves the volume of the bubble to nearly machine precision and accurately captures discontinuities in the pressure and modified conformation tensor across the interface, where there are jumps in density and viscosity.</div></div>","PeriodicalId":352,"journal":{"name":"Journal of Computational Physics","volume":"552 ","pages":"Article 114695"},"PeriodicalIF":3.8,"publicationDate":"2026-05-01","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"146025860","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}
京公网安备 11010802042870号
京ICP备2023020795号-1
扫码关注我们
求助内容：
应助结果提醒方式：