Robust Numerical Techniques Research Articles

Enhancing Accuracy and Efficiency in Stiff ODE Integration Using Variable Step Diagonal BBDF Approaches

Recent advancements in mathematical modelling have uncovered a growing number of systems exhibiting stiffness, a phenomenon that challenges the effectiveness of traditional numerical methods. Motivated by the need for more robust numerical techniques to address this issue, this paper presents an enhanced version of the Diagonally Block Backward Differentiation Formula (BBDF) that incorporates intermediate points, known as off-step points, to improve the accuracy and efficiency of solutions for stiff ordinary differential equations (ODEs). The new scheme leverages an adaptive step-size strategy to refine accuracy and efficiency between regular and off-grid integration steps. Theoretical analysis confirms that the proposed scheme is an A-stable and convergent method, as it satisfies the fundamental criteria of consistency, zero-stability, and A-stability. Numerical experiments on single and multivariable systems across varying time scales demonstrate significant improvements in solving stiff ODEs compared to existing techniques. Therefore, the new proposed method is an effective solver for stiff ODEs.

Phase-field approach for precise fracture tracking in anisotropic rocks: Integrating orthotropy-based energy decomposition and two-fold symmetric fracture toughness

Rock formations are known to exhibit material anisotropy, both in terms of elastic and fracture properties. This means that the fracture path in such formations is not a priori known but rather a complex unknown that requires robust numerical techniques to predict accurately. In this context, the phase-field model is considered particularly effective, provided that certain physical considerations are carefully adjusted to align with the physics of the problem. While addressing elastic anisotropy is well-established, the tension–compression asymmetry necessary to inhibit crack interpenetration in phase-field fracture models needs to account for the specific material anisotropy. Additionally, to accurately capture crack propagation, it is critical to simultaneously account for orientation-dependent fracture toughness in such materials. To address this, the present study employs an anisotropic phase-field model that integrates the generalized spectral decomposition proposed in the literature for orthotropic materials with a two-fold symmetric fracture toughness, to predict the fracture trajectories in rock-type samples under fixed mixed-mode loading ratios. While each of the two aspects has primarily been applied to model orthotropic plates under simple tensile and shearing loading conditions in the literature, here we study their applicability in complex loading scenarios. To this end, the experimental data from notched semi-circular specimens of Grimsel Granite undergoing complex mixed-mode loading obtained in our previous work is considered. We focus on two given mode-mixity ratios and perform numerical studies. Our results emphasize the importance of considering this generalized decomposition for phase-field modeling of fracturing in rock-type materials, particularly under loading conditions where the crack might otherwise be unrealistically driven into the compressive region. Although certain features are well captured by considering anisotropy in elasticity alone, our findings demonstrate that incorporating a two-fold symmetric fracture toughness proves to be advantageous for more precise tracking of the fracture path.

Physics-Inspired Quantum Simulation of Resonating Valence Bond States─A Prototypical Template for a Spin-Liquid Ground State.

Spin liquids─an emergent, exotic collective phase of matter─have garnered enormous attention in recent years. While experimentally many prospective candidates have been proposed and realized, theoretically modeling real materials that display such behavior may pose serious challenges due to the inherently high correlation content of such phases. Over the last few decades, the second-quantum revolution has been the harbinger of a novel computational paradigm capable of initiating a foundational evolution in computational physics. In this report, we strive to use the power of the latter to study a prototypical model, a spin-1/2-unit cell of a Kagome antiferromagnet. Extended lattices of such unit cells are known to possess a magnetically disordered spin-liquid ground state. We employ robust classical numerical techniques such as the density-matrix renormalization group (DMRG) to identify the nature of the ground state through a matrix-product state (MPS) formulation. We subsequently use the gained insight to construct an auxiliary Hamiltonian with reduced measurables and also design an ansatz that is modular and gate-efficient. With robust error-mitigation strategies, we are able to demonstrate that the said ansatz is capable of accurately representing the target ground state even on a real IBMQ backend within 1% accuracy in energy. Since the protocol is linearly scaling O(n) in the number of unit cells, gate requirements, and the number of measurements, it is straightforwardly extendable to larger Kagome lattices that can pave the way for efficient construction of spin-liquid ground states on a quantum device.

A fully coupled superficial runoff and soil erosion basin scale model with efficient time stepping

We present a numerical model of soil erosion at the basin scale that allows one to describe surface run-off without a priori identifying drainage zones, river beds and other water bodies. The model is based on robust semi-implicit numerical techniques and guarantees exact mass conservation and positivity of the surface and subsurface water layers. Furthermore, the method is equipped with a geostatistical preprocessor that can perform downscaling of data retrieved from digital databases at coarser resolutions. Numerical experiments on both idealized and realistic configurations demonstrate the effectiveness of the proposed method in reproducing transient high resolution features at a reduced computational cost and to reproduce correctly the main hydrographic features of the considered catchment. Furthermore, probabilistic forecasts can be carried out based on soil data maps automatically generated by the geostatistical preprocessor that are distributed among the available MPI ranks to carry out simulations independently thus reducing the total cost of the simulation. The numerical experiments show the ability of the model to provide robust estimates of water levels, discharge and of the order of magnitude of the total sediment yield.

Computational Irreversibility Analysis in Transient Flow of a Non-Newtonian Fluid through a Channel with Porous Walls

A second law analysis is explored to investigate the irreversibility properties in a transient, porous channel flow of a viscous, incompressible, and non-Newtonian fluid. The non-Newtonian fluid model is of a Generalized Newtonian Fluid type with no elastic properties but with shear-thinning viscosity. Additionally, given that the flow is non-isothermal, the viscosity is therefore expectedly also assumed to be temperature dependent. The porous channel is subjected to constant suction and injection of fluid through the walls. Computational solutions for the underlying fluid dynamical equations, based on robust finite difference numerical techniques, are developed and implemented in time and space. We demonstrate the effects of the embedded fluid flow and heat transfer parameters on the fluid velocity and temperature profiles. We also explore the competing effects of heat transfer irreversibility versus fluid friction irreversibility. The major observations are that, in the flow regions where the maximum velocity obtains, heat transfer irreversibility significantly dominates over the otherwise insignificant fluid friction irreversibility. It is also observed that, in those flow regions away from the region of maximum velocity, the opposite scenario obtains, and hence fluid friction irreversibility significantly dominates over heat transfer irreversibility. Along the channel walls, fluid friction irreversibility notably dominates over heat transfer irreversibility. The possibility that certain parameter choices may lead to reverse scenario is not discounted.

Physics-informed graph neural networks enhance scalability of variational nonequilibrium optimal control.

When a physical system is driven away from equilibrium, the statistical distribution of its dynamical trajectories informs many of its physical properties. Characterizing the nature of the distribution of dynamical observables, such as a current or entropy production rate, has become a central problem in nonequilibrium statistical mechanics. Asymptotically, for a broad class of observables, the distribution of a given observable satisfies a large deviation principle when the dynamics is Markovian, meaning that fluctuations can be characterized in the long-time limit by computing a scaled cumulant generating function. Calculating this function is not tractable analytically (nor often numerically) for complex, interacting systems, so the development of robust numerical techniques to carry out this computation is needed to probe the properties of nonequilibrium materials. Here, we describe an algorithm that recasts this task as an optimal control problem that can be solved variationally. We solve for optimal control forces using neural network ansatz that are tailored to the physical systems to which the forces are applied. We demonstrate that this approach leads to transferable and accurate solutions in two systems featuring large numbers of interacting particles.

Numerical daemons of hydrological models are summoned by extreme precipitation

Abstract. Hydrological models are usually systems of nonlinear differential equations for which no analytical solutions exist and thus rely on numerical solutions. While some studies have investigated the relationship between numerical method choice and model error, the extent to which extreme precipitation such as that observed during hurricanes Harvey and Katrina impacts numerical error of hydrological models is still unknown. This knowledge is relevant in light of climate change, where many regions will likely experience more intense precipitation. In this experiment, a large number of hydrographs are generated with the modular modeling framework FUSE (Framework for Understanding Structural Errors), using eight numerical techniques across a variety of forcing data sets. All constructed models are conceptual and lumped. Multiple model structures, parameter sets, and initial conditions are incorporated for generality. The computational cost and numerical error associated with each hydrograph were recorded. Numerical error is assessed via root mean square error and normalized root mean square error. It was found that the root mean square error usually increases with precipitation intensity and decreases with event duration. Some numerical methods constrain errors much more effectively than others, sometimes by many orders of magnitude. Of the tested numerical methods, a second-order adaptive explicit method is found to be the most efficient because it has both a small numerical error and a low computational cost. A small literature review indicates that many popular modeling codes use numerical techniques that were suggested by this experiment to be suboptimal. We conclude that relatively large numerical errors may be common in current models, highlighting the need for robust numerical techniques, in particular in the face of increasing precipitation extremes.

A Study of the COVID-19 Impacts on the Canadian Population

With the recent outbreak of COVID-19, the reach and scale of COVID-19 cases is top of mind for everyone and many research groups are actively monitoring and exploring the potential spread. A positive consequence of past epidemics and pandemics is that there are sound epidemiological compartmental modelling approaches that can effectively model disease spread. With minor changes to the underlying dynamical system of equations, many different strategies and situations can be explored. In particular, one such strategy of social distancing is top of mind for many Canadians as our political leaders, local businesses, and fellow Canadians promote and adopt this approach with the hopes that it will effectively `flatten the curve' and reduce or prevent further spread. In this paper, the baseline SIR model is introduced with its close counterpart, the SEIR model. Social distancing is modelled through the isolation of a subset of the susceptible population and comparative studies are performed considering a range in the proportion of individuals isolated. Robust and accurate numerical approximation techniques are used to simulate the pessimistic base case for which no preventative measures are taken and for various social distancing regimes. The results of social distancing are consolidated into two groups - those that flatten the curve and those that completely halt the disease spread. Mathematical formulations show that the turning point between these two regimes is when the effective reproductive rate, denoted $R_{e}$ , is equal to 1. Conclusions are made regarding the impacts and extent of the spread in relation to the severity of social distancing measures.

The space-time conservation element and solution element scheme for simulating two-phase flow in pipes

In this article, the space-time conservation element and solution element scheme is extended to simulate the unsteady compressible two-phase flow in pipes. The model is non-conservative and the governing equations consist of three equations, namely, two mass conservation equations for each phase and one mixture-momentum equation. In the third equation, the non-conservative source term appears, which describes the sum of gravitational and frictional forces. The presence of source term and two mass conservation equations in considered model offers difficulties in developing the accurate and robust numerical techniques. The suggested space-time conservation element and solution element numerical scheme resolves the volume-contact discontinuities efficiently. Furthermore, the modified central upwind scheme is also extended to solve the same two-phase flow model. The number of test problems is considered, and the results obtained by space-time conservation element and solution element scheme are compared with the solutions of modified central upwind scheme. The numerical results show better performance of the space-time conservation element and solution element method as compare to the modified central upwind scheme.

Numerical approximation of the 3D hydrostatic Navier–Stokes system with free surface

In this paper we propose a stable and robust strategy to approximate the 3D incompressible hydrostatic Euler and Navier–Stokes systems with free surface. Compared to shallow water approximation of the Navier–Stokes system, the idea is to use a Galerkin type approximation of the velocity field with piecewise constant basis functions in order to obtain an accurate description of the vertical profile of the horizontal velocity. Such a strategy has several advantages. It allowsto rewrite the Navier–Stokes equations under the form of a system of conservation laws with source terms,the easy handling of the free surface, which does not require moving meshes,the possibility to take advantage of robust and accurate numerical techniques developed in extensive amount for Shallow Water type systems.Compared to previous works of some of the authors, the three dimensional case is studied in this paper. We show that the model admits a kinetic interpretation including the vertical exchanges terms, and we use this result to formulate a robust finite volume scheme for its numerical approximation. All the aspects of the discrete scheme (fluxes, boundary conditions, ...) are completely described and the stability properties of the proposed numerical scheme (well-balancing, positivity of the water depth, ...) are discussed. We validate the model and the discrete scheme with some numerical academic examples (3D non stationary analytical solutions) and illustrate the capability of the discrete model to reproduce realistic tsunami waves propagation, tsunami runup and complex 3D hydrodynamics in a raceway.

Approximate Green's Function Coupled Cluster Method Employing Effective Dimension Reduction.

The Green's function coupled cluster (GFCC) method, originally proposed in the early 1990s, is a powerful many-body tool for computing and analyzing the electronic structure of molecular and periodic systems, especially when electrons of the system are strongly correlated. However, in order for the GFCC to become a method that may be routinely used in the electronic structure calculations, robust numerical techniques and approximations must be employed to reduce its extremely high computational overhead. In our recent studies, it has been demonstrated that the GFCC equations can be solved directly in the frequency domain using iterative linear solvers, which can be easily distributed in a massively parallel environment. In the present work, we demonstrate a successful application of model-order-reduction (MOR) techniques in the GFCC framework. Briefly speaking, for a frequency regime of interest that requires high-resolution descriptions of spectral function, instead of solving the GFCC linear equation of full dimension for every single frequency point of interest, an efficiently solvable linear system model of a reduced dimension may be built upon projecting the original GFCC linear system onto a subspace. From this reduced order model is obtained a reasonable approximation to the full dimensional GFCC linear equations in both interpolative and extrapolative spectral regions. Here, we show that the subspace can be properly constructed in an iterative manner from the auxiliary vectors of the GFCC linear equations at some selected frequencies within the spectral region of interest. During the iterations, the quality of the subspace, as well as the linear system model, can be systematically improved. The method is tested in this work in terms of the efficiency and accuracy of computing spectral functions for some typical molecular systems such as carbon monoxide, 1,3-butadiene, benzene, and adenine. To reach the same level of accuracy as that of the original GFCC method, the application of MOR in the GFCC method is able to significantly lower the original computational cost for the aforementioned molecules in designated frequency regimes. As a byproduct, the reduced order model obtained by this method is found to provide a high-quality initial guess, which improves the convergence rate for the existing iterative linear solver.

Unsteady advection-diffusion simulations in complex Discrete Fracture Networks with an optimization approach

It is widely recognized that the prediction of transport of contaminants in a fractured rock mass requires models that preserve several distinctive features of the inner fracture network, like heterogeneity and directionality; in this respect, Discrete Fracture Networks (DFNs) play a significant role. The solution of the associated equations would claim a high computational demand, that could be met only by using agile and robust numerical techniques. In this note a new numerical technique, fully validated from a mathematical standpoint, is applied to engineering problems, also introducing dispersion models for the description of non-stationary transport phenomena. The method results in a fast and scalable resolution tool, based on a PDE-constrained optimization approach designed to avoid mesh generation problems and allowing for transport simulations with an Eulerian approach. Examples are reported to show the quality of the solution obtained, even by using relatively coarse meshes and quite geometrically complex DFNs.

Differential Geometry in Edge Detection: Accurate Estimation of Position, Orientation and Curvature.

The vast majority of edge detection literature has aimed at improving edge recall and precision, with relatively few addressing the accuracy of edge orientation estimates which are often based on gradient. We show that first-order estimates of orientation can have significant error and this can be remedied by employing Third-Order estimates. This paper aims at estimating differential geometry attributes of an edge, namely, localization, orientation, and curvature, as well as edge topology, and develop robust numerical techniques in gray-scale and color images, applicable to a variety of popular edge detectors, such as gradient-based, gPb and SE. Second, a combinatorial model of edge grouping in a small neighborhood is developed to capture all geometrically consistent grouping called curvels, which establish: (i) edge topology in the form of potential links between an edge and other edges; (ii) an accurate curvature estimate for each possible grouping, whose performance is comparable to methods which use global and multi-scale methods; (iii) a more accurate localization of an edge. These have been evaluated using four distinct methodologies (i) traditional human annotated datasets; (ii) using coherence measure; (iii) stability analysis under visual perturbation, and (iv) utilitarian evaluation, and show meaningful improvements.

Quantifying credit portfolio losses under multi-factor models

ABSTRACTIn this work, we investigate the challenging problem of estimating credit risk measures of portfolios with exposure concentration under the multi-factor Gaussian and multi-factor t-copula models. It is well-known that Monte Carlo (MC) methods are highly demanding from the computational point of view in the aforementioned situations. We present efficient and robust numerical techniques based on the Haar wavelets theory for recovering the cumulative distribution function of the loss variable from its characteristic function. To the best of our knowledge, this is the first time that multi-factor t-copula models are considered outside the MC framework. The analysis of the approximation error and the results obtained in the numerical experiments section show a reliable and useful machinery for credit risk capital measurement purposes in line with Pillar II of the Basel Accords.

A conservative fully implicit algorithm for predicting slug flows

An accurate and predictive modelling of slug flows is required by many industries (e.g., oil and gas, nuclear engineering, chemical engineering) to prevent undesired events potentially leading to serious environmental accidents. For example, the hydrodynamic and terrain-induced slugging leads to unwanted unsteady flow conditions. This demands the development of fast and robust numerical techniques for predicting slug flows. The presented in this paper study proposes a multi-fluid model and its implementation method accounting for phase appearance and disappearance. The numerical modelling of phase appearance and disappearance presents a complex numerical challenge for all multi-component and multi-fluid models. Numerical challenges arise from the singular systems of equations when some phases are absent and from the solution discontinuity when some phases appear or disappear. This paper provides a flexible and robust solution to these issues. A fully implicit formulation described in this work enables to efficiently solve governing fluid flow equations. The proposed numerical method provides a modelling capability of phase appearance and disappearance processes, which is based on switching procedure between various sets of governing equations. These sets of equations are constructed using information about the number of phases present in the computational domain. The proposed scheme does not require an explicit truncation of solutions leading to a conservative scheme for mass and linear momentum. A transient two-fluid model is used to verify and validate the proposed algorithm for conditions of hydrodynamic and terrain-induced slug flow regimes. The developed modelling capabilities allow to predict all the major features of the experimental data, and are in a good quantitative agreement with them.

Heat Transfer in an L Shaped Porous Medium using FEM

Heat transfer in porous medium has been an intense field of study for many years that has tremendously improved our understanding about the heat flow behavior in porous region. Finite element method has been widely acknowledged to be one of the most robust numerical techniques that has proven its capability in handling the tedious and complex set of equations representing various phenomenon. The present research investigates the heat transfer inside a L shaped porous domain having heated from bottom surface and cooled at the top edge. Finite element method is used with the help of triangular element to divide the porous domain into smaller segments. The results are discussed with respect to various physical parameters affecting the heat transfer behaviour.

Numerical simulations of chaotic and complex spatiotemporal patterns in fractional reaction–diffusion systems

The generalized fractional reaction–diffusion equations which exist in the form of noninteger order partial differential equations have now found wide application for illustrating important and useful physical phenomena, such as subdiffusive and superdiffusive scenarios. The space fractional derivatives are defined in the Riesz sense on the intervals $$0<\alpha <1$$ and $$1<\alpha \le 2$$ . We propose robust numerical techniques, such as a spectral representation of the fractional Laplacian operator in conjunction with the exponential time differencing method, and present the equivalent relationship between the Riesz fractional derivative and fractional Laplacian operator. We apply these techniques to numerically solve a range of chaotic processes, such as the Chua’s equations, Rossler system, Lorenz and Lorenz-type systems. Simulation results revealed various complex and spatiotemporal chaos, spiral chaos, intermittent chaos and spots patterns in two-dimensional space fractional reaction–diffusion problems.

Quantifying Credit Portfolio Losses Under Multi-Factor Models

In this work, we investigate the challenging problem of estimating credit risk measures of portfolios with exposure concentration under the multi-factor Gaussian and multi-factor t-copula models. It is well-known that Monte Carlo (MC) methods are highly demanding from the computational point of view in the aforementioned situations. We present efficient and robust numerical techniques based on the Haar wavelets theory for recovering the cumulative distribution function (CDF) of the loss variable from its characteristic function. To the best of our knowledge, this is the first time that multi-factor t-copula models are considered outside the MC framework. The analysis of the approximation error and the results obtained in the numerical experiments section show a reliable and useful machinery for credit risk capital measurement purposes in line with Pillar II of the Basel Accords.

An immersed free-surface boundary treatment for seismic wave simulation

Finite-difference methods are popular for wave simulation within the seismic exploration community, thanks to their efficiency. However, difficulties arise when encountering complex topography due to the regular grid pattern of the finite-difference schemes. Despite alternatives that can handle the free surface with little effort, such as the spectral element or discontinuous Galerkin’s methods, incorporating a free-surface boundary condition within the finite-difference framework is still appealing, even at the cost of extra algorithm complexity and stronger requirement of computational resources. We present a free-surface boundary treatment within the finite-difference framework, belonging to the family of the immersed-boundary methods. Inherently, the presented boundary treatment is separated from the rest of the wave simulation, which makes it easy to be integrated in existing finite-difference codes. Specifically, we construct an extrapolation operator for each grid point above the free surface, if requested by the finite-difference stencil, to estimate its fictitious wavefield value at each time step. These operators are constructed only once and remain unchanged for all the time steps and source locations. The memory requirement of these operators is significant. Fortunately, grouping together multiple simulations concerning different source locations makes it possible to dilute the memory burden to a negligible level. Additionally, applying these operators incurs numerical noise, which may lead to long time instabilities. In such a scenario, additional numerical procedures, for instance, introducing artificial diffusion, are necessary to control the instability and obtain sensible simulation results. Successful applications of the presented boundary treatment to elastic-wave equations on domains with nontrivial topographies, in 2D and 3D, are presented. Robust and efficient numerical techniques to control high-frequency numerical noise remain to be investigated.

Physically based modeling in catchment hydrology at 50: Survey and outlook

Integrated, process-based numerical models in hydrology are rapidly evolving, spurred by novel theories in mathematical physics, advances in computational methods, insights from laboratory and field experiments, and the need to better understand and predict the potential impacts of population, land use, and climate change on our water resources. At the catchment scale, these simulation models are commonly based on conservation principles for surface and subsurface water flow and solute transport (e.g., the Richards, shallow water, and advection-dispersion equations), and they require robust numerical techniques for their resolution. Traditional (and still open) challenges in developing reliable and efficient models are associated with heterogeneity and variability in parameters and state variables; nonlinearities and scale effects in process dynamics; and complex or poorly known boundary conditions and initial system states. As catchment modeling enters a highly interdisciplinary era, new challenges arise from the need to maintain physical and numerical consistency in the description of multiple processes that interact over a range of scales and across different compartments of an overall system. This paper first gives an historical overview (past 50 years) of some of the key developments in physically based hydrological modeling, emphasizing how the interplay between theory, experiments, and modeling has contributed to advancing the state of the art. The second part of the paper examines some outstanding problems in integrated catchment modeling from the perspective of recent developments in mathematical and computational science.