2025-04-07 | | Total: 22
High-order discontinuous Galerkin spectral element methods (DGSEM) have received growing attention and development, especially in the regime of computational fluid dynamics in recent years. The inherent flexibility of the discontinuous Galerkin approach in handling non-conforming interfaces, such as those encountered in moving geometries or hp-refinement, presents a significant advantage for real-world simulations. Despite the well-established mathematical framework of DG methods, practical implementation challenges persist to boost performance and capability. Most previous studies only focus on certain choices of element shape or basis type in a structured mesh, although they have demonstrated the capability of DGSEM in complex flow simulations. This work discusses the low-cost and unified interface flux evaluation approaches for general spectral elements in unstructured meshes, alongside their implementations in the open-source spectral element framework, Nektar++. The initial motivation arises from the discretization of Helmholtz equations by the symmetric interior penalty method, in which the system matrix can easily become non-symmetric if the flux is not properly evaluated on non-conforming interfaces. We focus on the polynomial non-conforming case in this work but extending to the geometric non-conforming case is theoretically possible. Comparisons of different approaches, trade-offs, and performance of benchmark of our initial matrix-free implementation are also included, contributing to the broader discourse on high-performance spectral element method implementations.
This paper is concerned with inserting three-dimensional computer-aided design (CAD) geometries into meshes composed of hexahedral elements using a volume fraction representation. An adaptive procedure for doing so is presented. The procedure consists of two steps. The first step performs spatial acceleration using a k-d tree. The second step involves subdividing individual hexahedra in an adaptive mesh refinement (AMR)-like fashion and approximating the CAD geometry linearly (as a plane) at the finest subdivision. The procedure requires only two geometric queries from a CAD kernel: determining whether or not a queried spatial coordinate is inside or outside the CAD geometry and determining the closest point on the CAD geometry's surface from a given spatial coordinate. We prove that the procedure is second-order accurate for sufficiently smooth geometries and sufficiently refined background meshes. We demonstrate the expected order of accuracy is achieved with several verification tests and illustrate the procedure's effectiveness for several exemplar CAD geometries.
We establish Sobolev-Poincaré inequalities for piecewise W1,p functions over sequences of fairly general polytopic (thence also shape-regular simplicial and Cartesian) meshes in any dimension; amongst others, they cover the case of standard Poincaré inequalities for piecewise W1,p functions and can be useful in the analysis of nonconforming finite element discretizations of nonlinear problems. Crucial tools in their derivation are novel Sobolev-trace inequalities and Babu\v ska-Aziz inequalities with mixed boundary conditions. We provide estimates that are constant free, i.e., that are fully explicit with respect to the geometric properties of the domain and the underlying sequence of polytopic meshes.
Multiphysics and multiscale mathematical models enable the non-invasive study of cardiac function. These models often rely on simplifying assumptions that neglect certain biophysical processes to balance fidelity and computational cost. In this work, we propose an eikonal-based framework that incorporates mechano-calcium feedback -- the effect of mechanical deformation on calcium-troponin buffering -- while introducing only negligible computational overhead. To assess the impact of mechano-calcium feedback at the organ level, we develop a bidirectionally coupled cellular electromechanical model and integrate it into two cardiac multiscale frameworks: a monodomain-driven model that accounts for geometric feedback on electrophysiology and the proposed eikonal-based approach, which instead neglects geometric feedback. By ensuring consistent cellular model calibration across all scenarios, we isolate the role of mechano-calcium feedback and systematically compare its effects against models without it. Our results indicate that, under baseline conditions, mechano-calcium feedback has minimal influence on overall cardiac function. However, its effects become more pronounced in altered force generation scenarios, such as inotropic modulation. Furthermore, we demonstrate that the eikonal-based framework, despite omitting other types of mechano-electric feedback, effectively captures the role of mechano-calcium feedback at significantly lower computational costs than the monodomain-driven model, reinforcing its utility in computational cardiology.
In this paper, we extend our research concerning the standard and linearized monotonicity methods for the inverse problem of the time harmonic elastic wave equation and introduce the modification of these methods for noisy data. In more detail, the methods must provide consistent results when using noisy data in order to be able to perform simulations with real world data, e.g., laboratory data. We therefore consider the disturbed Neumann-to-Dirichlet operator and modify the bound of the eigenvalues in the monotonicity tests for reconstructing unknown inclusions with noisy data. In doing so, we show that there exists a noise level δ0 so that the inclusions are detected and their shape is reconstructed for all noise levels δ<δ0. Finally, we present some numerical simulations based on noisy data.
The paper presents a novel multimesh rational approximation scheme for the numerical solution of the (homogeneous) Dirichlet problem for the spectral fractional Laplacian. The scheme combines a rational approximation of the function λ↦λ−s with a set of finite element approximations of parameter-dependent non-fractional partial differential equations (PDEs). The key idea that underpins the proposed scheme is that each parametric PDE is numerically solved on an individually tailored finite element mesh. This is in contrast to the existing single-mesh approach, where the same finite element mesh is employed for solving all parametric PDEs. We develop an a posteriori error estimation strategy for the proposed rational approximation scheme and design an adaptive multimesh refinement algorithm. Numerical experiments show improvements in convergence rates compared to the rates for uniform mesh refinement and up to 10 times reduction in computational costs compared to the corresponding adaptive algorithm in the single-mesh setting.
In the present paper we consider the semiclassical magnetic Schrödinger equation, which describes the dynamics of charged particles under the influence of a electro-magnetic field. The solution of the time-dependent Schrödinger equation is approximated by a single Gaussian wave packet via the time-dependent Dirac--Frenkel variational principle. For the approximation we use ordinary differential equations of motion for the parameters of the variational solution and extend the second-order Boris algorithm for classical mechanics to the quantum mechanical case. In addition, we propose a modified version of the classical fourth order Runge--Kutta method. Numerical experiments explore parameter convergence and geometric properties. Moreover, we benchmark against the analytical solution of the Penning trap.
A space-discretization for the elastic flow of inextensible curves is devised and quasi-optimal convergence of the corresponding semi-discrete problem is proved for a suitable discretization of the nonlinear inextensibility constraint. Further a fully discrete time-stepping scheme that incorporates this constraint is proposed and unconditional stability and convergence of the discrete scheme are proved. Finally some numerical simulations are used to verify the obtained results experimentally.
We present a variationally consistent wrinkling model based on spectral decomposition of the stress tensor, providing a unified formulation that captures the three distinct membrane states. Compared to the previous strain-based spectral decomposition approach, the proposed model improves accuracy by satisfying the uniaxial tension condition from tension field theory and aligning with the mixed wrinkling criterion. It also demonstrates excellent performance under various loading conditions and offers enhanced generality by unifying strain-based, stress-based, and mixed criteria within a single framework. Beyond these improvements, the model retains the superior convergence properties of the previous approach, including the framework for the flexible inclusion or omission of residual compressive stiffness. This mitigates nonconvergence or singularities in slackening states. With these adjustments, new expressions for stress and constitutive tensors are consistently derived. Finally, extensive validation through analytical, numerical, and experimental benchmark tests highlights the robustness of the model. The results confirm its accuracy in capturing the mechanical response of wrinkled thin membranes, strong convergence properties, and value for advanced membrane wrinkling analysis.
Bayesian inverse problems are an important application for probabilistic solvers of partial differential equations: when fully resolving numerical error is computationally infeasible, probabilistic solvers can be used to consistently model the error and propagate it to the posterior. In this work, the performance of the random mesh finite element method (RM-FEM) is investigated in a Bayesian inverse setting. We show how interpolation error negatively affects the RM-FEM posterior, and how these negative effects can be diminished. In scenarios where FEM is biased for a quantity of interest, we find that RM-FEM struggles to accurately model this bias.
We analyze a first-order exponential wave integrator (EWI) for the nonlinear Schrödinger equation (NLSE) with a singular potential locally in L2, which might be locally unbounded. The typical example is the inverse power potential such as the Coulomb potential, which is the most fundamental potential in quantum physics and chemistry. We prove that, under the assumption of L2-potential and H2-initial data, the L2-norm convergence of the EWI is, roughly, first-order in one dimension (1D) and two dimensions (2D), and 34-order in three dimensions (3D). In addition, under a stronger integrability assumption of Lp-potential for some p>2 in 3D, the L2-norm convergence increases to almost 34+3(12−1p) order if p≤125 and becomes first-order if p>125. In particular, our results show, to the best of our knowledge for the first time, that first-order L2-norm convergence can be achieved when solving the NLSE with the Coulomb potential in 3D. The key advancements are the use of discrete (in time) Strichartz estimates, which allow us to handle the loss of integrability due to the singular potential that does not belong to L∞, and the more favorable local truncation error of the EWI, which requires no (spatial) smoothness of the potential. Extensive numerical results in 1D, 2D, and 3D are reported to confirm our error estimates and to show the sharpness of our assumptions on the regularity of the singular potentials.
Multirate integration is an increasingly relevant tool that enables scientists to simulate multiphysics systems. Existing multirate methods are designed for equations whose fast and slow variables can be linearly separated using additive or component-wise partitions. However, in realistic applications, this assumption is not always valid. Building on the recently developed class of nonlinearly partitioned Runge-Kutta (NPRK) methods, we develop a framework for multirate NPRK (MR-NPRK) that allows for arbitrary nonlinear splittings of the evolution operator. We discuss order conditions, formalize different types of coupling between timescales, and analyze joint linear stability of MR-NPRK methods. We then introduce a class of 2nd- and 3rd-order methods, referred to as ``implicitly-wrapped'' multirate methods, that combine a user-specified explicit method for integrating the fast timescale with several slow implicit stages. These methods are designed to be algorithmically simple with low memory costs and minimal operator evaluations. Lastly, we conduct numerical experiments to validate our proposed methods and show the benefits of multirating a nonlinear partition.
We study a system of diffusing point particles in which any triplet of particles reacts and is removed from the system when the relative proximity of the constituent particles satisfies a predefined condition. Proximity-based reaction conditions of this kind are commonly used in particle-based simulations of chemical kinetics to mimic bimolecular reactions, those involving just two reactants, and have been extensively studied. The rate at which particles react within the system is determined by the reaction condition and particulate diffusion. In the bimolecular case, analytic relations exist between the reaction rate and the distance at which particles react allowing modellers to tune the rate of the reaction within their simulations by simply altering the reaction condition. However, generalising proximity-based reaction conditions to trimolecular reactions, those involving three particles, is more complicated because it requires understanding the distribution of the closest diffusing particle to a point in the vicinity of a spatially dependent absorbing boundary condition. We find that in this case the evolution of the system is described by a nonlinear partial integro-differential equation with no known analytic solution, which makes it difficult to relate the reaction rate to the reaction condition. To resolve this, we use singular perturbation theory to obtain a leading-order solution and show how to derive an approximate expression for the reaction rate. We then use finite element methods to quantify the higher-order corrections to this solution and the reaction rate, which are difficult to obtain analytically. Leveraging the insights gathered from this analysis, we demonstrate how to correct for the errors that arise from adopting the approximate expression for the reaction rate, enabling for the construction of more accurate particle-based simulations than previously possible.
This study investigates enhancing option pricing by extending the Black-Scholes model to include stochastic volatility and interest rate variability within the Partial Differential Equation (PDE). The PDE is solved using the finite difference method. The extended Black-Scholes model and a machine learning-based LSTM model are developed and evaluated for pricing Google stock options. Both models were backtested using historical market data. While the LSTM model exhibited higher predictive accuracy, the finite difference method demonstrated superior computational efficiency. This work provides insights into model performance under varying market conditions and emphasizes the potential of hybrid approaches for robust financial modeling.
Surface acoustic wave (SAW) devices are widely used in modern communication equipment and SAW equations describe the critical physical processes of acoustic-electric conversion in SAW devices. It is very challenging to numerically solve such equations, since they are typically three dimensional problems defined on unbounded domains. In this paper, we first use the perfectly matched layer method to truncate the unbounded domain and then propose a finite element tearing and interconnecting algorithm for the truncated equations based on the periodic structure of the truncated domain. We also design an effective solver for the ill-conditioned linear system of the Lagrange multipliers arising from discretization. Several numerical results are performed to demonstrate the efficiency of the proposed algorithm.
We develop efficient and high-order accurate solvers for the Helmholtz equation on complex geometry. The schemes are based on the WaveHoltz algorithm which computes solutions of the Helmholtz equation by time-filtering solutions of the wave equation. The approach avoids the need to invert an indefinite matrix which can cause convergence difficulties for many iterative solvers for indefinite Helmholtz problems. Complex geometry is treated with overset grids which use Cartesian grids throughout most of the domain together with curvilinear grids near boundaries. The basic WaveHoltz fixed-point iteration is accelerated using GMRES and also by a deflation technique using a set of precomputed eigenmodes. The solution of the wave equation is solved efficiently with implicit time-stepping using as few as five time-steps per period, independent of the mesh size. The time-domain solver is adjusted to remove dispersion errors in time and this enables the use of such large time-steps without degrading the accuracy. When multigrid is used to solve the implicit time-stepping equations, the cost of the resulting WaveHoltz scheme scales linearly with the total number of grid points N (at fixed frequency) and is thus optimal in CPU-time and memory usage as the mesh is refined. A simple rule-of-thumb formula is provided to estimate the number of points-per-wavelength required for a p-th order accurate scheme which accounts for pollution (dispersion) errors. Numerical results are given for problems in two and three space dimensions, to second and fourth-order accuracy, and they show the potential of the approach to solve a wide range of large-scale problems.
This paper investigates the semi-discrete optimal transport (OT) problem with entropic regularization. We characterize the solution using a governing, well-posed ordinary differential equation (ODE). This naturally yields an algorithm to solve the problem numerically, which we prove has desirable properties, notably including global strong convexity of a value function whose Hessian must be inverted in the numerical scheme. Extensive numerical experiments are conducted to validate our approach. We compare the solutions obtained using the ODE method with those derived from Newton's method. Our results demonstrate that the proposed algorithm is competitive for problems involving the squared Euclidean distance and exhibits superior performance when applied to various powers of the Euclidean distance. Finally, we note that the ODE approach yields an estimate on the rate of convergence of the solution as the regularization parameter vanishes, for a generic cost function.
MaRDI Open Interfaces is a software project aimed at improving reuse and interoperability in Scientific Computing by alleviating the difficulties of crossing boundaries between different programming languages, in which numerical packages are usually implemented, and of switching between multiple implementations of the same mathematical problem. The software consists of a set of formal interface specifications for common Scientific Computing tasks, as well as a set of loosely coupled libraries that facilitate implementing these interfaces or adapting existing implementations for multiple programming languages and handle data marshalling automatically without sacrificing performance, enabling users to use different implementations without significant code efforts. The software has high reuse potential due to aim to solve general numerical problems.
We introduce and compute solutions of a many-body model of the electronic properties of twisted bilayer graphene which systematically accounts for the effects of structural relaxation. We model mechanical relaxation by coupling linear elasticity to a stacking energy that penalizes disregistry. Minimizers of the resulting functional are then input into a tight-binding model of twisted bilayer graphene, from which a single-particle continuum moiré-scale (Bistritzer-MacDonald-like) model is systematically derived. We then project this model together with a Coulomb electron-electron interaction term into the single-particle model's flat moiré bands. We numerically compute Hartree-Fock ground states of this model, comparing the relative energies of competing many-body ground states.
Phakopsora pachyrhizi is a biotrophic fungus that requires living plant tissue to survive in the wild for extended periods. The fungus causes Asian rust and costs billions of US dollars every year for its control in South American soy production. Despite the regulatory measure that prohibits the cultivation of soybeans in some months of the year (soybean-free period) in Brazil, the presence of soybean production areas in neighboring countries, such as Paraguay and Bolivia, can help the survival of the pathogen between crop seasons. It is known that P. pachyrhizi urediniospores can be spread/transported thousands of kilometres by the wind. In this context, the objective of this work was to develop a mathematical model to simulate the atmospheric transport of P. pachyrhizi urediniospores from the west, Paraguay/Bolivia to Paraná Brazil, through storms coming from cold fronts. The transport of urediniospores was modeled by a diffusive-convective-reactive equation. Wind direction and the velocity of cold fronts that crossed Paraná state between October 2018 and February 2019 were used. For validation, real cases of rust occurrence in the Paraná state informed by the Anti-rust Consortium Portal in the season 2018/19 were used. The results confirm mathematically that urediniospores from infected areas located in a country on the west can be transported and deposited on the east, in the state of Paraná. The first case of soybean rust in Paraná state/Brazil was registered 10 days after the first cold front, suggesting that the transported and deposited urediniospores were still viable for host infection. This work reinforces the importance of the establishment of the soybean-free period in other soybean-producing countries. It will also provide a better understanding of the fungus dispersion system, potentially enabling the correct use of fungicides.
This paper studies properties of fixed points of generalised Extra-gradient (GEG) algorithms applied to min-max problems. We discuss connections between saddle points of the objective function of the min-max problem and GEG fixed points. We show that, under appropriate step-size selections, the set of saddle points (Nash equilibria) is a subset of stable fixed points of GEG. Convergence properties of the GEG algorithm are obtained through a stability analysis of a discrete-time dynamical system. The results and benefits when compared to existing methods are illustrated through numerical examples.
Unconstrained optimization problems become more common in scientific computing and engineering applications with the rapid development of artificial intelligence, and numerical methods for solving them more quickly and efficiently have been getting more attention and research. Moreover, an efficient method to minimize all kinds of objective functions is urgently needed, especially the nonsmooth objective function. Therefore, in the current paper, we focus on proposing a novel numerical method tailored for unconstrained optimization problems whether the objective function is smooth or not. To be specific, based on the variational procedure to refine the gradient and Hessian matrix approximations, an efficient quadratic model with 2n constrained conditions is established. Moreover, to improve the computational efficiency, a simplified model with 2 constrained conditions is also proposed, where the gradient and Hessian matrix can be explicitly updated, and the corresponding boundedness of the remaining 2n−2 constrained conditions is derived. On the other hand, the novel numerical method is summarized, and approximation results on derivative information are also analyzed and shown. Numerical experiments involving smooth, derivative blasting, and non-smooth problems are tested, demonstrating its feasibility and efficiency. Compared with existing methods, our proposed method can efficiently solve smooth and non-smooth unconstrained optimization problems for the first time, and it is very easy to program the code, indicating that our proposed method not also has great application prospects, but is also very meaningful to explore practical complex engineering and scientific problems.