This paper proposes a quasi-asymptotic-preserving hybrid discontinuous Galerkin (HDG-QAP) scheme for the resolution of highly anisotropic diffusion problems. The HDG-QAP scheme introduces an auxiliary unknown that serves to capture the information on the dominant diffusion scale. We show that it is well-posed for any ε > 0, with ε being a small constant and 1/ ε characterizing the anisotropy strength, and that its solution is bounded uniformly in ε . At this point, the standard HDG procedure, namely, the static condensation on the numerical trace, turns out to violate these uniform bounds. Instead, we show that the static condensation on the auxiliary unknown does lead to similar bounds uniform in ε , but its resolution is costly in terms of computational effort. Therefore, we propose a relaxation method, called the HDG-QAP Uzawa iteration, to overcome this challenge, in which each iteration is fast to compute. We show that the HDG-QAP Uzawa iteration converges for any ε > 0, but also for ε = 0. Finally, we provide some numerical examples to confirm the findings, and in particular to show that the proposed HDG-QAP Uzawa iteration works well even with a severe anisotropy ε = 10 −15 where the quality of the numerical solutions is unaffected by the anisotropy strength.
The invariant energy quadratization (IEQ) framework has emerged as a powerful tool for designing energy-stable numerical schemes for gradient flow problems. However, for the Allen–Cahn equation, it remains an open problem whether the IEQ-type schemes can rigorously preserve the discrete maximum bound principle (MBP)—a critical property ensuring the physical admissibility of solutions. While numerical evidence has consistently suggested that IEQ-based schemes tend to preserve the MBP in practice, no theoretical justification has been established to date. In this work, we develop and analyze a fully discrete, second-order accurate numerical scheme based on the IEQ framework that rigorously satisfies both the discrete energy stability and the MBP. The key innovation lies in a new treatment of the auxiliary variable: instead of using second-order extrapolation, we employ a numerically bounded approximation derived from a first-order MBP-preserving scheme, ensuring that the nonlinear coefficient remains controlled. We prove that the resulting scheme is energy stable, and provably satisfies a discrete MBP under a mild and explicitly characterized time-step condition. In addition, we establish optimal-order error estimates for the fully discrete system. Extensive 2D and 3D numerical experiments are conducted to validate the theoretical analysis and demonstrate the accuracy, efficiency, and structure-preserving properties of the proposed method. To the best of our knowledge, this is the first theoretical result showing that an IEQ-based scheme with second-order time discretization can rigorously preserve the MBP, thereby addressing a fundamental open question in the field.
In this paper, a stable and arbitrary high order spectral volume (SV) method is proposed for solving linear hyperbolic conservation laws on unstructured quadrilateral meshes. The SV scheme is constructed with subdivision points being zeros of a class of parameterized polynomials. A unified proof for the $L^2$ stability of the SV method is provided under conditions that the parameter $c>-\frac{1}{k(k+1)}$ and the underlying mesh is an $h^{1+\gamma}, \gamma\ge 1$ parallelogram mesh. The optimal a priori error estimate of SV methods is established. Numerical experiments are presented to verify all theoretical findings.
Abstract. Froth flotation in a column is a widely used unit operation in mineral processing, wastewater treatment, and other applications. The flotation process selectively separates finely divided hydrophobic materials (valuable minerals or ores; repelled by water) from hydrophilic (slimes or gangue; attracted to water), where both are suspended in a viscous fluid. A flotation column roughly functions as follows: gas is introduced close to the bottom and generates bubbles that rise through the continuously injected pulp that contains the solid particles. The hydrophobic particles attach to the bubbles, forming foam or froth (the concentrate) that is removed through a launder. The hydrophilic particles do not attach to bubbles, but normally settle to the bottom, and are removed continuously. Additional wash water, injected close to the top, may assist with the rejection of entrained impurities and increase froth stability. A recently formulated partial differential equation model [R. Bürger, S. Diehl, M.C. Martí, Y. Vásquez, IMA J. Appl. Math. 87 (2022) 1151–1190] describes the process by a pair of degenerate parabolic PDEs with discontinuous flux for the volume fractions of bubbles and hydrophilic solid particles as functions of height and time. An extension of that model is presented, which includes the effect of wash water to be injected into the froth as well as the transport of an arbitrary number of components (such as slimes or chemical reagents) with the liquid. The numerical scheme for bubble and particle volume fractions is extended to simulate percentages representing the liquid components, which are proven to remain nonnegative and sum up to one. In addition, a theory of desired steady states of the flotation column is outlined. It is proven that the condition of “positive bias” in a determined zone of the flotation column (i.e., net downward flow of water) coincides with the mathematically derived condition for the existence of a stationary bubble concentration profile, including a stable froth layer. It is demonstrated how steady-state solutions to the governing model can be constructed and conditions for their existence can be conveniently mapped through so-called “operating charts.” Numerical simulations are presented.
We construct a non-polynomial local discontinuous Galerkin (LDG) scheme for the prescribed mean curvature equation to approximate boundary gradient blow-up solutions and obtain error estimates.
The solution to the elastodynamic equations with mixed boundary conditions exhibits a singular behavior where the boundary conditions change. We obtain detailed expansions of the singularities of the Dirichlet trace of the solution and the traction on the boundary. The results imply quasi-optimal estimates for piecewise polynomial approximations. The results are applied to hp and graded versions of the time domain boundary element method. Numerical examples illustrate the theoretical results in 2d. They confirm the expected quasi-optimal convergence rates and the singular behavior of the numerical approximations.
We present a two-dimensional free-boundary Hele-Shaw model to describe certain aspects of eukaryotic cell motility on substrates. The key ingredients of this model are Darcy's law for the over-damping motion of the cytoplasm as a confined viscous droplet, with a particular boundary condition for the Young-Laplace equation. This particular condition describes the active force induced by the cytoskeleton, which can be generated in cells either by actin polymerization against the membrane or by actomyosin contraction of cortical filaments, which adhere to the membrane. In addition, in this model we take into account the friction exerted by the membrane on the substrate. First, we study the linear stability of the steady state and prove that above a threshold, the disk is linearly unstable. This analysis highlights the stabilizing effect of undercooling. Then, using a bifurcation argument, we prove the existence of traveling waves that describe a persistent motion in cell migration and justify the relevance of the model. We emphasize in particular that the friction exerted by the membrane on the substrate has a stabilizing effect on cell dynamics.
We are investigating the numerical solution to the 2D time-harmonic Maxwell equations in the presence of a classical medium and a metamaterial, that is with sign-changing coefficients. As soon as the problem has a (unique) solution, we are able to build a converging numerical approximation based on the finite element method, for which there is no constraint on the meshes related to the sign-changing behavior. To that aim, we use Lagrange finite elements to approximate the scalar potentials appearing in the Helmholtz decomposition of the vector-valued electromagnetic fields. Convergence in strong norm is proven for the fields. Numerical examples illustrate the theory.
This article analyzes hybrid high-order method for space discretization and backward Euler and Crank-Nicolson schemes for time discretization of the nonlinear extended Fisher-Kolmogorov and the Fisher-Kolmogorov equations. The critical parameter gamma > 0 is incorporated in the stabilization term of the HHO method. Error estimate of order & Oscr;(h(k+1) + Delta t) (resp. & Oscr;(h(k+1) + (Delta t)(2)) in the energy-norm for the backward Euler (resp. Crank-Nicolson) scheme is obtained when polynomials of order k+2 (resp. k) with k >= 0 are utilized to approximate the exact solution in the interior of the polygon and its traces on the boundary of the polygon (resp. normal derivative on the mesh faces). The HHO discretization for the Fisher-Kolmogorov equation, that is, when the parameter gamma = 0, leads to a convergence rate of & Oscr;(h(k+2)) in the space variable. The results of the numerical experiments validate the theoretical results.
In this work, we derive a hyperbolic system of dispersive equations for the numerical simulation of coastal waves with improved dispersive properties and admitting an exact energy conservation equation. This system is derived with the assumption of a moderate non-linearity and of a correction coefficient close to 1. This system contains the same non-linear terms as the Serre-Green-Naghdi equations, which are obtained in the limit where the Mach number tends to zero. The assumptions are only used to neglect non-linear terms related to the improvement of dispersive properties. The bathymetry can be included with a mild-slope hypothesis. On this basis, we propose an energy-stable numerical scheme relying on a splitting between the hyperbolic and dispersive parts of the model. The stability of the method is achieved through the discrete dissipation of the energy balance specific to each step. We also establish the existence of soliton solutions for this model. Numerical simulations are proposed to highlight the dispersive properties of the model, as well as the dissipative character of the scheme.
This paper initiates a mathematical investigation of a PDE model for the transport of high voltage direct current via overhead lines. We prove the existence of infinitely many solutions, give necessary conditions for existence, explicitly compute the continuum of all radial solutions, and develop a new numerical algorithm for this problem.
We study finite-difference approximations of the Poisson-Boltzmann (PB) electrostatic energy functional of ionic concentrations and electric displacements constrained by Gauss' law and the ionic mass conservation, and a class of local algorithms for minimizing the finite-difference discretized such energy functional. We prove that the discrete Boltzmann distributions characterize the finite-difference minimizer and obtain the uniform bounds and optimal error estimates in maximum norm for such a minimizer. The local algorithm is an iteration over all the grid boxes that locally minimizes the energy by updating the concentrations and displacement one grid box at a time, keeping Gauss' law and the mass conservation satisfied. A new local algorithm with a shift is constructed for minimizing the Poisson electrostatic energy (the part of the PB energy without ionic concentrations) with a variable dielectric coefficient. We prove the convergence of these local algorithms and present numerical tests to demonstrate the results of our analysis.
In this work, we develop an efficient numerical method for solving 3D Maxwell's equations in non-cylindrical coaxial cables. The main challenge arises from the elongated geometry of the computational domain, which induces strong anisotropy between the longitudinal direction (along the cable) and the transverse directions (within the cross-sections). This leads to the use of highly anisotropic meshes, where the longitudinal mesh size is much larger than the transverse one. Our objective is to design a numerical scheme that is explicit in the longitudinal direction, with a CFL stability condition depending only on the longitudinal mesh size. In a previous work, we achieved this for cylindrical cables by employing prismatic edge elements, 1D quadrature for longitudinal mass lumping, and a hybrid explicit/implicit time discretization. The present paper extends this approach to non-cylindrical cables, addressing several new difficulties with the following key ingredients: (1) representing the cable as a deformation of a reference cylindrical cable and employing mapping techniques between the physical and reference domains; (2) using an anisotropic space discretization that combines an interior penalty discontinuous Galerkin (IPDG) method in the transverse directions with a conforming finite element method in the longitudinal direction; (3) utilizing prismatic edge elements on a prismatic mesh of the reference cable; and (4) adapting the construction of the hybrid explicit-implicit time discretization to the new structure of the semi-discrete problem. From a theoretical perspective, the main difficulty lies in the stability analysis, which requires extending and adapting standard techniques for DG methods in space and energy methods in time.
In this work, we focus on developing a new class of numerical methods able to handle both quasineutrality and charge separation in plasmas. At large temporal and spatial scales, plasmas tend to be quasineutral, meaning that the local net charge density is nearly zero. However, when the scale at which one observes the plasma dynamics is smaller than the characteristic distance over which the electric field and charges are typically screened, then quasineutrality breaks down. In such regimes, standard numerical methods face severe stability constraints, rendering them practically unusable. To address this issue, in this work, we introduce and analyze a new class of finite volume penalized-IMEX Runge-Kutta methods for the Euler-Poisson system, specifically designed to handle the quasineutral limit. We show that, these proposed schemes are uniformly stable with respect to the Debye length and degenerate into high order methods as the quasineutral limit is approached. Several numerical tests confirm that this new class of methods exhibits the desired properties.
We introduce and analyze a new mixed finite element method for the stationary model arising from the coupling of the Brinkman-Forchheimer and Darcy equations. While the original unknowns are given by the velocities and pressures of the more and less permeable porous media, our approach is based on the introduction of the Brinkman-Forchheimer pseudostress as a further variable, which allows us to eliminate the respective pressure from the formulation. Nevertheless, this latter unknown, along with other variables of physical interest, such as the velocity gradient, vorticity, and the stress tensor, can be accurately recovered afterwards by means of postprocessing formulae that depend mainly on the pseudostress, all of which constitutes one of the most distinctive feature of the present strategy. Next, aiming to perform a proper treatment of the transmission conditions, the traces on the interface, of both the Brinkman Forchheimer velocity and the Darcy pressure, are also incorporated as auxiliary unknowns. Thus, the resulting fully-mixed variational formulation can be seen as a nonlinear perturbation of, in turn, a twofold perturbed saddle point operator equation. Additionally, the diagonal feature of some of the bilinear forms involved, facilitates the proof of their corresponding inf-sup conditions. Then, the fixed-point strategy arising from a linearization of the Forchheimer term, along with suitable abstract results exploiting the aforementioned structure and the classical Banach theorem, are employed to prove the existence and uniqueness of a solution under a suitable small-data assumption, both for the fully-mixed variational formulation and for the discrete scheme arising from the associated Galerkin system. In particular, Raviart Thomas and piecewise polynomial subspaces of the lowest degree for the domain unknowns, as well as continuous piecewise linear polynomials for the interface ones, constitute a feasible choice. Under this selection of spaces, momentum is conserved in both the Brinkman-Forchheimer and Darcy equations whenever the external forces belong to the piecewise constants, thus yielding another relevant characteristic of our approach. Optimal error estimates and associated rates of convergence are established. Finally, several numerical results illustrating the good performance of the method and confirming the theoretical findings, are reported.
In this work, we propose an improved discretization, in terms of stability and accuracy, for the incompressible two-phase Darcy flows in a heterogeneous porous medium with discontinuous capillary forces. For this purpose, the total velocity formulation of the model is used. The coupled system is composed of a degenerate parabolic equation for the non-wetting phase and a pressure equation for the total velocity. We combine a positive Vertex Approximation Gradient (VAG) type scheme for the gradient fluxes with a hybrid upwinding of the mobilities. This approach entails a maximum principle on the saturations, which remain in their physical ranges. Energy estimates are obtained by selecting key approximations of the fluxes. These stability results allow to prove the existence of discrete solutions. Numerical experiments on complex test-cases show the robustness of the new approach in terms of the accuracy as well as the nonlinear convergence. Comparison to the usual phase potential upwinding approach and to a previous hybrid upwinding scheme are also provided.
In this paper, we propose and analyze a strongly mass-conservative numerical scheme for the coupled Navier-Stokes and Darcy-Forchheimer system in both two and three spatial dimensions. The two subproblems are coupled through physically relevant interface conditions, including mass conservation, balance of normal forces, and the Beavers-Joseph-Saffman condition. We employ a staggered discontinuous Galerkin method for the Navier-Stokes equations and use standard mixed finite elements for the Darcy-Forchheimer problem. The proposed formulation incorporates the interface conditions directly, without introducing Lagrange multipliers on the interface or artificial numerical fluxes on the mesh skeleton. As a consequence, although discontinuous Galerkin elements are used in the free-flow region, the resulting discrete velocity field is globally H(div)-conforming across the entire domain. In particular, the incompressibility constraint is satisfied exactly in the free-flow region, thereby yielding strong mass conservation over the entire computational domain. Under a suitable small-data assumption, we establish the well-posedness of the resulting nonlinear discrete system. Owing to the exact preservation of mass conservation, the proposed scheme exhibits a pressure-robust behavior, in the sense that the velocity approximation is insensitive to pressure effects. Numerical experiments are presented to illustrate the stability and robustness of the method, including its performance in regimes involving small viscosity, large pressure, and limited solution regularity.
We revisit the higher-order Localized Orthogonal Decomposition variant by Maier [SIAM J. Numer. Anal. 59 (2021) 1067-1089] based on nonconforming constraints (discontinuous finite element spaces) and introduce a new variant based on conforming constraints (continuous finite elements), putting both approaches in a general unified framework. We propose a new localization strategy that is suitable for both approaches and offers a new perspective on the localization of LOD in general. We fully analyze the strategy for linear scalar elliptic problems and discuss extensions to the Helmholtz equation and the Gross-Pitaevskii eigenvalue problem. Numerical examples are presented that provide valuable comparisons between conforming and nonconforming constraints.
In this paper, we are concerned with two-scale integrators for the non-relativistic Klein-Gordon (NRKG) equation with a dimensionless parameter 0 < e << 1, which is inversely proportional to the speed of light. The highly oscillatory property in time of this model corresponds to the parameter e and the equation in the form of partial derivative(tt)u-Delta/epsilon(2)u + 1/epsilon(4)u + lambda/epsilon(2)f(u) = 0 has a factor 1/epsilon(2) in front of the nonlinearity which means that this part becomes strong when e is small. These two aspects bring significantly numerical burdens in designing numerical methods. We propose a class of two-scale integrators which is constructed based on some reformulations to the system, Fourier pseudo-spectral method and exponential integrators. Two practical integrators up to order three and four are constructed by using some symmetric conditions and the stiff order conditions of implicit exponential integrators. The convergence of the obtained integrators is rigorously studied, and it is shown that the uniform accuracy in time is O(h(3)) and O(h(4)) for the time stepsize h. The near energy conservation over long times is also established for the multi-stage integrators by using modulated Fourier expansions. Numerical results on a NRKG equation show that the proposed integrators have high accuracy, excellent long time energy conservation and competitive efficiency.
We conduct a comprehensive convergence and error analysis of the second- and third-order unconditionally energy stable convex splitting Runge-Kutta (CSRK) methods for H−1 gradient flows with typical forms of free energy. Through the energy structure inherent to gradient flows, we are able to derive uniform-in-time bounds of the numerical solution in the H1, H2, and L6 norms. In turn, these functional bounds enable us to derive the associated estimates for the nonlinear error terms. Meanwhile, motivated by the fact that the diffusion coefficients are diagonally dominated in the CSRK numerical systems, the convergence results become available, based on a stage-by-stage analysis for the error evolutionary equations. The Cahn–Hilliard and phase-field crystal equations are two examples in the theoretical analysis. We also numerically compute some convergence results to validate the theorems proposed in this paper. This work deepens the theoretical foundation of CSRK methods and provides robust analytical tools for their application to conserved gradient flows.
