Enhanced geothermal systems (EGS) involve strongly coupled, advection-dominated flow and heat transfer in fractured porous media. Conventional models typically assume local thermal equilibrium with a single effective fluid temperature or, at best, an averaged pore-fluid temperature, so the thermal evolution of injected cold fluid is only inferred indirectly. In this work, we develop a local thermal non-equilibrium (LTNE) model that explicitly resolves the temperature of injected fluid as it moves through the reservoir and exchanges heat with the hot rock and resident fluid. The key ingredient is a concentration variable that tracks the injected fluid and induces a three-way LTNE coupling among rock, resident-fluid, and injected-fluid temperatures. This framework distinguishes, at the continuum scale, how newly injected fluid parcels are heated by conductive and convective exchange, and predicts production-well temperatures without relying on bulk averages. To discretize the resulting nonlinear, advection-dominated system, we employ an enriched Galerkin (EG) finite element method for Darcy flow, temperature, and concentration, providing local mass conservation with relatively few degrees of freedom. We further design a flux-corrected transport (FCT) strategy for the EG concentration and temperature equations to enforce a discrete maximum principle and suppress nonphysical oscillations while preserving local conservation. Time integration uses an IMPES-type splitting combined with a strong-stability-preserving Runge–Kutta scheme. Numerical experiments for fractured EGS problems show that the proposed LTNE–EG–FCT framework captures injected-fluid heating paths and thermal breakthrough behavior not resolved by standard single-temperature or averaged LTNE models.

We propose a new kind of localized shock capturing for continuous (CG) and discontinuous Galerkin (DG) discretizations of hyperbolic conservation laws. The underlying framework of dissipation-based weighted essentially nonoscillatory (WENO) stabilization for high-order CG and DG approximations was introduced in our previous work. In this general framework, Hermite WENO (HWENO) reconstructions are used to calculate local smoothness sensors that determine the appropriate amount of artificial viscosity for each cell. In the original version, candidate polynomials for WENO averaging are constructed using the derivative data from von Neumann neighbors. We upgrade this standard 'cell-cell' reconstruction procedure by using WENO polynomials associated with mesh vertices as candidate polynomials for cell-based WENO averaging. The Hermite data of individual cells is sent to vertices of those cells, after which vertex-averaged HWENO data is sent back to cells containing the vertices. The new 'cell-vertex' averaging procedure includes the data of vertex neighbors without explicitly adding them to the reconstruction stencils. It mitigates mesh imprinting and can also be used in classical HWENO limiters for DG methods. The second main novelty of the proposed approach is a quadrature-driven distribution of artificial viscosity within high-order finite elements. Replacing the linear quadrature weights by their nonlinear WENO-type counterparts, we concentrate shock-capturing dissipation near discontinuities while minimizing it in smooth portions of troubled cells. This redistribution of WENO stabilization preserves the total dissipation rate within each cell and improves local shock resolution without relying on subcell decomposition techniques. Numerical experiments in one and two dimensions demonstrate substantial improvements in accuracy and robustness for high-order elements, providing a compelling alternative to standard cell-cell reconstructions and subcell shock-capturing schemes.

We present an invariant-domain-preserving (IDP) treatment of nonconforming interfaces for Legendre–Gauss–Lobatto Discontinuous Galerkin Spectral Element Methods (LGL-DGSEM) with adaptive mesh refinement (AMR) on Cartesian meshes. The proposed methodology extends recently developed convex limiting and graph-viscosity frameworks for DGSEM to meshes containing hanging nodes. Starting from a conservative mortar formulation, we derive low-order interface fluxes that satisfy the requirements of invariant-domain-preserving discretizations. To avoid the excessive diffusion associated with fully connected mortar couplings, a sparsification strategy based on LGL subcell characteristic functions is introduced, yielding compact interface stencils. The resulting mortar fluxes remain conservative, reduce to the standard conforming formulation on matching interfaces, and naturally fit into graph-viscosity-based low-order schemes used for convex limiting. The proposed construction provides the missing ingredient required to combine high-order DGSEM discretizations, invariant-domain-preserving limiting, and adaptive mesh refinement within a unified framework for nonlinear hyperbolic conservation laws. We provide numerical verifications of the properties of the proposed scheme and run challenging simulations that require positivity limiting and shock-capturing.

Accurate prediction of shallow water flows relies on precise bottom topography data, yet direct bathymetric surveys are expensive and time-consuming. In contrast, remote sensing platforms such as radar or satellite altimetry provide accurate free surface observations. This disparity motivates a data-driven reconstruction strategy: invert the shallow water equations to estimate the bathymetry that yields the best fit to the governing dynamics. We introduce a new direct reconstruction technique that extracts bathymetric features from widely available free surface measurements. The underlying inverse problem of determining an unknown bathymetry profile from observed wave elevations is inherently ill-posed. Small perturbations in the data may lead to large deviations in the reconstructed topography, and discontinuities or sharp gradients further exacerbate instability. To stabilize the inversion, we formulate an optimal-control problem, wherein a cost functional penalizes deviations between simulated and measured free surface elevation while enforcing a state equation for the flow dynamics. To suppress noise and preserve sharp depth variations, the framework is augmented with L^1 regularization and total variation denoising. These sparsity-promoting terms encourage piecewise-smooth solutions, allowing changes in the bathymetry to be captured without excessive smoothing. Numerical experiments on synthetic noisy data and discontinuous bathymetry demonstrate robust performance in reconstructing unknown bathymetry.

We discretize the M_1 model of radiative transfer using continuous finite elements and propose a tailor-made monolithic convex limiting (MCL) procedure for enforcing physical realizability. The M_1 system of nonlinear balance laws for the zeroth and first moments of a probability distribution function is derived from the linear Boltzmann equation and equipped with an entropy-based closure for the second moment. To ensure hyperbolicity and physical admissibility, evolving moments must stay in an invariant domain representing a convex set of realizable states. We first construct a low-order method that is provably invariant domain preserving (IDP). Introducing intermediate states that represent spatially averaged exact solutions of homogeneous Riemann problems, we prove that these so-called bar states are realizable in any number of space dimensions. This key auxiliary result enables us to show the IDP property of a fully discrete scheme with a diagonally implicit treatment of reactive terms. To achieve high resolution, we add nonlinear correction terms that are constrained using a two-step MCL algorithm. In the first limiting step, local bounds are imposed on each conserved variable to avoid spurious oscillations and maintain positivity of the scalar-valued zeroth moment (particle density). The second limiting step constrains the magnitude of the vector-valued first moment to be realizable. The flux-corrected finite element scheme is provably IDP. Its ability to prevent nonphysical behavior while attaining high-order accuracy in smooth regions is verified in a series of numerical tests. The developed methodology provides a robust simulation tool for dose calculation in radiotherapy.

Simulating infiltration in porous media using Richards' equation remains computationally challenging due to its parabolic structure and nonlinear coefficients. While a wide range of numerical methods for differential equations have been applied over the past several decades, basic higher-order numerical methods often fail to preserve physical bounds on water pressure and saturation, leading to spurious oscillations and poor iterative solver convergence. Instead, low-order, bound-preserving methods have been preferred. The combination of mass lumping and relative permeability upwinding preserves bounds but degrades accuracy to first order in space. Flux-corrected transport is a high-resolution numerical technique designed for combining the bound-preserving property of low-order schemes with the accuracy of high-order methods, by blending the two methods through limited anti-diffusive fluxes. In this work, we extend flux-corrected transport schemes to the nonlinear, degenerate parabolic structure of Richards' equation, verify attainment of second-order convergence on unstructured meshes, and demonstrate applications to stormwater management infrastructure.

We review some recent advances in the field of element-based algebraic stabilization for continuous finite element discretizations of nonlinear hyperbolic problems. The main focus is on multidimensional convex limiting techniques designed to constrain antidiffusive element contributions rather than fluxes. We show that the resulting schemes can be interpreted as residual distribution methods. Two kinds of convex limiting can be used to enforce the validity of generalized discrete maximum principles in this context. The first approach has the structure of a localized flux-corrected transport (FCT) algorithm, in which the computation of a low-order predictor is followed by an antidiffusive correction stage. The second option is the use of a monolithic convex limiting (MCL) procedure at the level of spatial semi-discretization. In both cases, inequality constraints are imposed on scalar functions of intermediate states that are required to stay in convex invariant sets.

We equip a high-order continuous Galerkin discretization of a general hyperbolic problem with a nonlinear stabilization term and introduce a new methodology for enforcing preservation of invariant domains. The amount of shock-capturing artificial viscosity is determined by a smoothness sensor that measures deviations from a weighted essentially nonoscillatory (WENO) reconstruction. Since this kind of dissipative stabilization does not guarantee that the nodal states of the finite element approximation stay in a convex admissible set, we adaptively constrain deviations of these states from intermediate cell averages. The representation of our scheme in terms of such cell averages makes it possible to apply convex limiting techniques originally designed for positivity-preserving discontinuous Galerkin (DG) methods. Adapting these techniques to the continuous Galerkin setting and using Bernstein polynomials as local basis functions, we prove the invariant domain preservation property under a time step restriction that can be significantly weakened by using a flux limiter for the auxiliary cell averages. The close relationship to DG-WENO schemes is exploited and discussed. All algorithmic steps can be implemented in a matrix-free and hardware-aware manner. The effectiveness of the new element-based limiting strategy is illustrated by numerical examples.

We introduce a Lagrangian nodal discontinuous Galerkin (DG) hydrodynamics method for solving multi-dimensional hyperbolic systems. By incorporating an adaptation of Zalesak’s flux-corrected transport algorithm, we combine a first-order positivity-preserving scheme with a higher-order target discretization. This results in a flux-corrected Lagrangian DG scheme that ensures both global positivity preservation and second-order accuracy for the element averages of specific volume. The correction factors for flux limiting are derived from specific volume and applied to all components of the solution vector. We algebraically evolve the volumes of mesh elements using a discrete version of the geometric conservation law (GCL). The application of a limiter to the GCL fluxes is equivalent to moving the mesh using limited nodal velocities. Additionally, we equip our method with a locally bound-preserving slope limiter to effectively suppress spurious oscillations. Nodal velocity and external forces are computed using a multidirectional approximate Riemann solver to maintain conservation of momentum and total energy in vertex neighborhoods. Using linear finite elements together with a second-order time integrator ensures that the GCL is satisfied at the same order of accuracy. The results for standard test problems demonstrate the stability and superb shock-capturing capabilities of our scheme.

We present a deterministic framework for proton therapy dose calculation based on finite element discretizations of the energy-dependent M_1 moment model. The nonlinear M_1 system is derived from the Fokker–Planck equation for charged particles and closed using an entropy-based approximation of the second moment. Energy is treated as a pseudo-time coordinate. The zeroth and first moments of the proton fluence are evolved backward in energy. To ensure hyperbolicity and physical admissibility, we employ a monolithic convex limiting (MCL) strategy. Representing the standard continuous Galerkin discretization in terms of auxiliary `bar' states, we construct a nonlinear scheme that is provably invariant domain preserving (IDP) w.r.t. convex realizable sets consisting of all admissible states. The realizability of the bar states is enforced using the MCL technology for homogeneous hyperbolic systems. The forcing induced by stiff scattering is incorporated using Strang-type operator splitting. We use an explicit strong-stability-preserving Runge–Kutta method for the radiation transport subproblem and exact integration in the forcing steps, which guarantees the IDP property. The deposited dose is defined as the integral of a weighted zeroth moment over a bounded energy range. It is accumulated during the backward-in-energy evolution. Numerical experiments demonstrate that the proposed Strang-MCL method produces accurate and physically consistent dose distributions.

Dense particle suspensions are promising heat transfer fluids for next-generation Concentrated Solar Power (CSP) receivers, enabling operating temperatures above 800 degrees C. However, accurate modeling of the rheological behavior of granular flows is essential for reliable computational fluid dynamics (CFD) simulations. In this study, we develop and assess numerical methodologies for simulating dense suspensions pertinent to CSP applications. Our computational framework is based on Direct Numerical Simulation (DNS), augmented by lubrication force models to resolve detailed particle-particle and particle-wall interactions at volume fractions exceeding 50%. We conducted a systematic series of simulations across a range of volume fractions to establish a robust reference dataset. Validation was performed via a numerical viscometer configuration, permitting direct comparison with theoretical predictions and established benchmark results. Subsequently, the viscometer arrangement was generalized to a periodic cubic domain, serving as a representative volume element for CSP systems. Within this framework, effective viscosities were quantified independently through wall force measurements and energy dissipation fitting. The close agreement between these two approaches substantiates the reliability of the results. Based on these findings, effective viscosity tables were constructed and fitted using polynomial and piecewise-smooth approximations. These high-accuracy closure relations are suitable for incorporation into large-scale, non-Newtonian CFD models for CSP plant design.

We propose a way to maintain strong consistency and perform error analysis in the context of dissipation-based WENO stabilization for continuous and discontinuous Galerkin discretizations of conservation laws. Following Kuzmin and Vedral (J. Comput. Phys. 487:112153, 2023) and Vedral (arXiv preprint arXiv:2309.12019), we use WENO shock detectors to determine appropriate amounts of low-order artificial viscosity. In contrast to existing WENO methods, our approach blends candidate polynomials using residual-based nonlinear weights. The shock-capturing terms of our stabilized Galerkin methods vanish if residuals do. This enables us to achieve improved accuracy compared to weakly consistent alternatives. As we show in the context of steady convection-diffusion-reaction (CDR) equations, nonlinear local projection stabilization terms can be included in a way that preserves the coercivity of local bilinear forms. For the corresponding Galerkin-WENO discretization of a CDR problem, we rigorously derive a priori error estimates. Additionally, we demonstrate the stability and accuracy of the proposed method through one- and two-dimensional numerical experiments for hyperbolic conservation laws and systems thereof. The numerical results for representative test problems are superior to those obtained with traditional WENO schemes, particularly in scenarios involving shocks and steep gradients.

We introduce an unfitted Nitsche finite element method with a new ghost-penalty stabilization based on local projection of the solution gradient. The proposed ghost-penalty operator is straightforward to implement, ensures algebraic stability, provides an implicit extension of the solution beyond the physical domain, and stabilizes the numerical method for problems dominated by transport phenomena. This paper presents both a sharp interface version of the method and an alternative diffuse interface formulation designed to avoid integration over implicitly defined embedded surfaces. A complete numerical analysis of the sharp interface version is provided. The results of several numerical experiments support the theoretical analysis and illustrate the performance of both variants of the method.

In this work, we use the monolithic convex limiting (MCL) methodology to enforce relevant inequality constraints in implicit finite element discretizations of the compressible Euler equations. In this context, preservation of invariant domains follows from positivity preservation for intermediate states of the density and internal energy. To avoid spurious oscillations, we additionally impose local maximum principles on intermediate states of the density, velocity components, and specific total energy. For the backward Euler time stepping, we show the invariant domain preserving (IDP) property of the fully discrete MCL scheme by constructing a fixed-point iteration that meets the requirements of a Krasnoselskii-type theorem. Our iterative solver for the nonlinear discrete problem employs a more efficient fixed-point iteration. The matrix of the associated linear system is a robust low-order Jacobian approximation that exploits the homogeneity property of the flux function. The limited antidiffusive terms are treated explicitly. We use positivity preservation as a stopping criterion for nonlinear iterations. The first iteration yields the solution of a linearized semi-implicit problem. This solution possesses the discrete conservation property but is generally not IDP. Further iterations are performed if any non-IDP states are detected. The existence of an IDP limit is guaranteed by our analysis. To facilitate convergence to steady-state solutions, we perform adaptive explicit underrelaxation at the end of each time step. The calculation of appropriate relaxation factors is based on an approximate minimization of nodal entropy residuals. The performance of proposed algorithms and alternative solution strategies is illustrated by the convergence history for standard two-dimensional test problems.

We show that finite element discretizations of incompressible flow problems can be designed to ensure preservation/dissipation of kinetic energy not only globally but also locally. In the context of equal-order (piecewise-linear) interpolations, we prove the validity of a semi-discrete energy inequality for a quadraturebased approximation to the nonlinear convective term, which we combine with the Becker-Hansbo pressure stabilization. An analogy with entropy-stable algebraic flux correction schemes for the compressible Euler equations and the shallow water equations yields a weak 'bounded variation' estimate from which we deduce the semi-discrete Lax-Wendroff consistency and convergence towards dissipative weak solutions. The results of our numerical experiments for standard test problems confirm that the method under investigation is non-oscillatory and exhibits optimal convergence behavior.

We discretize the Vlasov-Poisson system using conservative semi-Lagrangian (CSL) discontinuous Galerkin (DG) schemes that are asymptotic preserving (AP) in the quasi-neutral limit. The proposed method (CSLDG) relies on two key ingredients: the CSLDG discretization and a reformulated Poisson equation (RPE). The use of the CSL formulation ensures local mass conservation and circumvents the Courant-Friedrichs-Lewy condition, while the DG method provides high-order accuracy for capturing fine-scale phase space structures of the distribution function. The RPE is derived by the Poisson equation coupled with moments of the Vlasov equation. The synergy between the CSLDG and RPE components makes it possible to obtain reliable numerical solutions, even when the spatial and temporal resolution might not fully resolve the Debye length. We rigorously prove that the proposed method is asymptotically stable, consistent and satisfies AP properties. Moreover, its efficiency is maintained across non-quasi-neutral and quasi-neutral regimes. These properties of our approach are essential for accurate and robust numerical simulation of complex electrostatic plasmas. Several numerical experiments verify the accuracy, stability and efficiency of the proposed CSLDG schemes.

We introduce a new multimesh finite element method for direct numerical simulation of incompressible particulate flows. The proposed approach falls into the category of overlapping domain decomposition / Chimera / overset grid meshes. In addition to calculating the velocity and pressure of the fictitious fluid on a fixed background mesh, we solve the incompressible Navier-Stokes equations on body-fitted submeshes that are attached to moving particles. The submesh velocity and pressure are used to calculate the hydrodynamic forces and torques acting on the particles. The coupling with the background velocity and pressure is enforced via (i) Robin-type boundary conditions for an Arbitrary-Lagrangian-Eulerian (ALE) formulation of the submesh problems and (ii) a Dirichlet-type distributed interior penalty term in the weak form of the background mesh problem. The implementation of the weak Dirichlet-Robin coupling is discussed in the context of discrete projection methods and finite element discretizations. Detailed numerical studies are performed for standard test problems involving fixed and moving immersed objects. A comparison of Chimera results with those produced by fictitious boundary methods illustrates significant gains in the accuracy of drag and lift approximations.

In this paper, we develop monolithic limiting techniques for enforcing nonlinear stability constraints in enriched Galerkin (EG) discretizations of nonlinear scalar hyperbolic equations. To achieve local mass conservation and gain control over the cell averages, the space of continuous (multi-)linear finite element approximations is enriched with piecewise-constant functions. The resulting spatial semi-discretization has the structure of a variational multiscale method. For linear advection equations, it is inherently stable but generally not bound preserving. To satisfy discrete maximum principles and ensure entropy stability in the nonlinear case, we use limiters adapted to the structure of our locally conservative EG method. The cell averages are constrained using a flux limiter, while the nodal values of the continuous component are constrained using a clip-and-scale limiting strategy for antidiffusive element contributions. The design and analysis of our new algorithms build on recent advances in the fields of convex limiting and algebraic entropy fixes for finite element methods. In addition to proving the claimed properties of the proposed approach, we conduct numerical studies for two-dimensional nonlinear hyperbolic problems. The numerical results demonstrate the ability of our limiters to prevent violations of the imposed constraints, while preserving the optimal order of accuracy in experiments with smooth solutions.

We investigate the consistency and convergence of flux-corrected finite element approximations in the context of nonlinear hyperbolic conservation laws. In particular, we focus on a monolithic convex limiting approach and prove a Lax–Wendroff-type theorem for the corresponding semi-discrete problem. A key component of our analysis is the use of a weak estimate on bounded variation, which follows from the semi-discrete entropy stability property of the method under investigation. For the Euler equations of gas dynamics, we prove the weak convergence of the flux-corrected finite element scheme to a dissipative weak solution under the assumption that a gas stays in physically reasonable region, i.e., the density is uniformly bounded away from zero, and the energy is uniformly bounded from above. Furthermore, if a strong solution exists, the sequence of numerical approximations converges strongly to the strong solution.

We propose a combination of machine learning and flux limiting for property-preserving subgrid scale modeling in the context of flux-limited finite volume methods for the one-dimensional shallow-water equations. The numerical fluxes of a conservative target scheme are fitted to the coarse-mesh averages of a monotone fine-grid discretization using a neural network to parametrize the subgrid scale components. To ensure positivity preservation and the validity of local maximum principles, we use a flux limiter that constrains the intermediate states of an equivalent fluctuation form to stay in a convex admissible set. The results of our numerical studies confirm that the proposed combination of machine learning with monolithic convex limiting produces meaningful closures even in scenarios for which the network was not trained.