A numerical methodology for enforcing maximum principles and the non-negative constraint for transient diffusion equations

Transient diffusion equations arise in many branches of engineering and applied sciences (e.g., heat transfer and mass transfer), and are parabolic partial differential equations. It is well-known that, under certain assumptions on the input data, these equations satisfy important mathematical properties like maximum principles and the non-negative constraint, which have implications in mathematical modeling. However, existing numerical formulations for these types of equations do not, in general, satisfy maximum principles and the non-negative constraint. In this paper, we present a methodology for enforcing maximum principles and the non-negative constraint for transient anisotropic diffusion equation. The method of horizontal lines (also known as the Rothe method) is applied in which the time is discretized first. This results in solving steady anisotropic diffusion equation with decay equation at every discrete time level. The proposed methodology for transient anisotropic diffusion equation will satisfy maximum principles and the non-negative constraint on general computational grids, and with no additional restrictions on the time step. We illustrate the performance and accuracy of the proposed formulation using representative numerical examples. We also perform numerical convergence of the proposed methodology. For comparison, we also present the results from the standard single-field semi-discrete formulation and the results from a popular software package, which all will violate maximum principles and the non-negative constraint

A scalable variational inequality approach for flow through porous media models with pressure-dependent viscosity

Mathematical models for flow through porous media typically enjoy the so-called maximum principles, which place bounds on the pressure field. It is highly desirable to preserve these bounds on the pressure field in predictive numerical simulations, that is, one needs to satisfy discrete maximum principles (DMP). Unfortunately, many of the existing formulations for flow through porous media models do not satisfy DMP. This paper presents a robust, scalable numerical formulation based on variational inequalities (VI), to model non-linear flows through heterogeneous, anisotropic porous media without violating DMP. VI is an optimization technique that places bounds on the numerical solutions of partial differential equations. To crystallize the ideas, a modification to Darcy equations by taking into account pressure-dependent viscosity will be discretized using the lowest-order Raviart-Thomas (RT0) and Variational Multi-scale (VMS) finite element formulations. It will be shown that these formulations violate DMP, and, in fact, these violations increase with an increase in anisotropy. It will be shown that the proposed VI-based formulation provides a viable route to enforce DMP. Moreover, it will be shown that the proposed formulation is scalable, and can work with any numerical discretization and weak form. Parallel scalability on modern computational platforms will be illustrated through strong-scaling studies, which will prove the efficiency of the proposed formulation in a parallel setting. Algorithmic scalability as the problem size is scaled up will be demonstrated through novel static-scaling studies. The performed static-scaling studies can serve as a guide for users to be able to select an appropriate discretization for a given problem size

On mesh restrictions to satisfy comparison principles, maximum principles, and the non-negative constraint: Recent developments and new results

This paper concerns with mesh restrictions that are needed to satisfy several important mathematical properties -- maximum principles, comparison principles, and the non-negative constraint -- for a general linear second-order elliptic partial differential equation. We critically review some recent developments in the field of discrete maximum principles, derive new results, and discuss some possible future research directions in this area. In particular, we derive restrictions for a three-node triangular (T3) element and a four-node quadrilateral (Q4) element to satisfy comparison principles, maximum principles, and the non-negative constraint under the standard single-field Galerkin formulation. Analysis is restricted to uniformly elliptic linear differential operators in divergence form with Dirichlet boundary conditions specified on the entire boundary of the domain. Various versions of maximum principles and comparison principles are discussed in both continuous and discrete settings. In the literature, it is well-known that an acute-angled triangle is sufficient to satisfy the discrete weak maximum principle for pure isotropic diffusion. An iterative algorithm is developed to construct simplicial meshes that preserves discrete maximum principles using existing open source mesh generators. Various numerical examples based on different types of triangulations are presented to show the pros and cons of placing restrictions on a computational mesh. We also quantify local and global mass conservation errors using representative numerical examples, and illustrate the performance of metric-based meshes with respect to mass conservation

A monolithic multi-time-step computational framework for first-order transient systems with disparate scales

Developing robust simulation tools for problems involving multiple mathematical scales has been a subject of great interest in computational mathematics and engineering. A desirable feature to have in a numerical formulation for multiscale transient problems is to be able to employ different time-steps (multi-time-step coupling), and different time integrators and different numerical formulations (mixed methods) in different regions of the computational domain. We present two new monolithic multi-time-step mixed coupling methods for first-order transient systems. We shall employ unsteady advection-diffusion-reaction equation with linear decay as the model problem, which offers several unique challenges in terms of non-self-adjoint spatial operator and rich features in the solutions. We shall employ the dual Schur domain decomposition technique to handle the decomposition of domain into subdomains. Two different methods of enforcing compatibility along the subdomain interface will be used in the time discrete setting. A systematic theoretical analysis (which includes numerical stability, influence of perturbations, bounds on drift along the subdomain interface) will be performed. The first coupling method ensures that there is no drift along the subdomain interface but does not facilitate explicit/implicit coupling. The second coupling method allows explicit/implicit coupling with controlled (but non-zero) drift in the solution along the subdomain interface. Several canonical problems will be solved to numerically verify the theoretical predictions, and to illustrate the overall performance of the proposed coupling methods. Finally, we shall illustrate the robustness of the proposed coupling methods using a multi-time-step transient simulation of a fast bimolecular advective-diffusive-reactive system

A deep learning modeling framework to capture mixing patterns in reactive-transport systems

Prediction and control of chemical mixing are vital for many scientific areas such as subsurface reactive transport, climate modeling, combustion, epidemiology, and pharmacology. Due to the complex nature of mixing in heterogeneous and anisotropic media, the mathematical models related to this phenomenon are not analytically tractable. Numerical simulations often provide a viable route to predict chemical mixing accurately. However, contemporary modeling approaches for mixing cannot utilize available spatial-temporal data to improve the accuracy of the future prediction and can be compute-intensive, especially when the spatial domain is large and for long-term temporal predictions. To address this knowledge gap, we will present in this paper a deep-learning (DL) modeling framework applied to predict the progress of chemical mixing under fast bimolecular reactions. This framework uses convolutional neural networks (CNN) for capturing spatial patterns and long short-term memory (LSTM) networks for forecasting temporal variations in mixing. By careful design of the framework -- placement of non-negative constraint on the weights of the CNN and the selection of activation function, the framework ensures non-negativity of the chemical species at all spatial points and for all times. Our DL-based framework is fast, accurate, and requires minimal data for training

Physics-Informed Machine Learning Models for Predicting the Progress of Reactive-Mixing

This paper presents a physics-informed machine learning (ML) framework to construct reduced-order models (ROMs) for reactive-transport quantities of interest (QoIs) based on high-fidelity numerical simulations. QoIs include species decay, product yield, and degree of mixing. The ROMs for QoIs are applied to quantify and understand how the chemical species evolve over time. First, high-resolution datasets for constructing ROMs are generated by solving anisotropic reaction-diffusion equations using a non-negative finite element formulation for different input parameters. Non-negative finite element formulation ensures that the species concentration is non-negative (which is needed for computing QoIs) on coarse computational grids even under high anisotropy. The reactive-mixing model input parameters are a time-scale associated with flipping of velocity, a spatial-scale controlling small/large vortex structures of velocity, a perturbation parameter of the vortex-based velocity, anisotropic dispersion strength/contrast, and molecular diffusion. Second, random forests, F-test, and mutual information criterion are used to evaluate the importance of model inputs/features with respect to QoIs. Third, Support Vector Machines (SVM) and Support Vector Regression (SVR) are used to construct ROMs based on the model inputs. Then, SVR-ROMs are used to predict scaling of QoIs. Qualitatively, SVR-ROMs are able to describe the trends observed in the scaling law associated with QoIs. Fourth, the scaling law's exponent dependence on model inputs/features are evaluated using $k$-means clustering. Finally, in terms of the computational cost, the proposed SVM-ROMs and SVR-ROMs are $\mathcal{O}(10^7)$ times faster than running a high-fidelity numerical simulation for evaluating QoIs

Modeling flow in porous media with double porosity/permeability: A stabilized mixed formulation, error analysis, and numerical solutions

The flow of incompressible fluids through porous media plays a crucial role in many technological applications such as enhanced oil recovery and geological carbon-dioxide sequestration. The flow within numerous natural and synthetic porous materials that contain multiple scales of pores cannot be adequately described by the classical Darcy equations. It is for this reason that mathematical models for fluid flow in media with multiple scales of pores have been proposed in the literature. However, these models are analytically intractable for realistic problems. In this paper, a stabilized mixed four-field finite element formulation is presented to study the flow of an incompressible fluid in porous media exhibiting double porosity/permeability. The stabilization terms and the stabilization parameters are derived in a mathematically and thermodynamically consistent manner, and the computationally convenient equal-order interpolation of all the field variables is shown to be stable. A systematic error analysis is performed on the resulting stabilized weak formulation. Representative problems, patch tests and numerical convergence analyses are performed to illustrate the performance and convergence behavior of the proposed mixed formulation in the discrete setting. The accuracy of numerical solutions is assessed using the mathematical properties satisfied by the solutions of this double porosity/permeability model. Moreover, it is shown that the proposed framework can perform well under transient conditions and that it can capture well-known instabilities such as viscous fingering

Modelling fluid deformable surfaces with an emphasis on biological interfaces

Fluid deformable surfaces are ubiquitous in cell and tissue biology, including lipid bilayers, the actomyosin cortex, or epithelial cell sheets. These interfaces exhibit a complex interplay between elasticity, low Reynolds number interfacial hydrodynamics, chemistry, and geometry, and govern important biological processes such as cellular traffic, division, migration, or tissue morphogenesis. To address the modelling challenges posed by this class of problems, in which interfacial phenomena tightly interact with the shape and dynamics of the surface, we develop a general continuum mechanics and computational framework for fluid deformable surfaces. The dual solid-fluid nature of fluid deformable surfaces challenges classical Lagrangian or Eulerian descriptions of deforming bodies. Here, we extend the notion of Arbitrarily Lagrangian-Eulerian (ALE) formulations, well-established for bulk media, to deforming surfaces. To systematically develop models for fluid deformable surfaces, which consistently treat all couplings between fields and geometry, we follow a nonlinear Onsager formalism according to which the dynamics minimize a Rayleighian functional where dissipation, power input and energy release rate compete. Finally, we propose new computational methods, which build on Onsager's formalism and our ALE formulation, to deal with the resulting stiff system of higher-order of partial differential equations. We apply our theoretical and computational methodology to classical models for lipid bilayers and the cell cortex. The methods developed here allow us to formulate/simulate these models for the first time in their full three-dimensional generality, accounting for finite curvatures and finite shape changes.Comment: 58 pages, 21 figures, submitted for publication to the Journal of Fluid Mechanic

Inference, prediction and optimization of non-pharmaceutical interventions using compartment models: the PyRoss library

PyRoss is an open-source Python library that offers an integrated platform for inference, prediction and optimisation of NPIs in age- and contact-structured epidemiological compartment models. This report outlines the rationale and functionality of the PyRoss library, with various illustrations and examples focusing on well-mixed, age-structured populations. The PyRoss library supports arbitrary structured models formulated stochastically (as master equations) or deterministically (as ODEs) and allows mid-run transitioning from one to the other. By supporting additional compartmental subdivision ad libitum, PyRoss can emulate time-since-infection models and allows medical stages such as hospitalization or quarantine to be modelled and forecast. The PyRoss library enables fitting to epidemiological data, as available, using Bayesian parameter inference, so that competing models can be weighed by their evidence. PyRoss allows fully Bayesian forecasts of the impact of idealized NPIs by convolving uncertainties arising from epidemiological data, model choice, parameters, and intrinsic stochasticity. Algorithms to optimize time-dependent NPI scenarios against user-defined cost functions are included. PyRoss's current age-structured compartment framework for well-mixed populations will in future reports be extended to include compartments structured by location, occupation, use of travel networks and other attributes relevant to assessing disease spread and the impact of NPIs. We argue that such compartment models, by allowing social data of arbitrary granularity to be combined with Bayesian parameter estimation for poorly-known disease variables, could enable more powerful and robust prediction than other approaches to detailed epidemic modelling. We invite others to use the PyRoss library for research to address today's COVID-19 crisis, and to plan for future pandemics.Comment: Code and updates at https://github.com/rajeshrinet/pyross 75 pages, 14 figures, and 1 tabl

Computational Methods for Multi-Scale Temporal Problems: Algorithms, Analysis, and Numerical Experiments

A major challenge in numerical simulation of most natural phenomena is the presence of disparate temporal and spatial scales. Capturing all the fine features can be computationally prohibitive. Hence, development of efficient and accurate multi-scale numerical algorithms has gained immense attention from engineers and scientists. Typically, a single numerical method cannot efficiently capture all the aforementioned features. Due to the assumptions made in construction of numerical methods and mathematical models, the range of applicability to various length and time-scales is often limited. A direction in resolving this issue is to apply different numerical methods in different regions of the computational domain. This strategy enables computation of necessary details as desired by the user. In this work, we propose numerical methodologies based on domain partitioning techniques that allow different time-steps and time-integrators in different regions of the computational domain. The first problem of interest is elastodynamics, which can pose various temporal scales in impact, contact and wave propagation problems. A monolithic (strong) coupling algorithm based on non-overlapping domain partitioning is proposed. The proposed algorithm is based on the theory of differential/algebraic equations and its numerical stability, energy conservation and accuracy is studied in detail. Following these findings, we extend this algorithm to advection-diffusion-reaction problems. The proposed algorithm proves useful especially in cases where the relative strength of the involved processes changes dramatically with respect to spatial coordinates. Numerical stability and accuracy of this method is studied and its application to fast bimolecular chemical reactions is showcased. Further on, we confine our attention to single and multiple-relaxation-time lattice Boltzmann methods for the advection-diffusion equation and study their performance in preserving the maximum principle and the non-negative constraint. Finally, a computational framework based on overlapping domain decomposition methods is proposed. This framework is designed for advection-diffusion problems and allows coupling of the finite element method and lattice Boltzmann methods with different time-steps and grid sizes. Additionally, a new method for enforcing the Dirichlet and Neumann boundary conditions on the numerical solution from the lattice Boltzmann method is proposed. This method is based on maximization of entropy and ensures non-negativity of the discrete distributions on the boundary of the domain. We study the performance of this framework through numerical experiments and showcase its application to fast and equilibrium chemical reactions.Civil and Environmental Engineering, Department o