Standard Discretization Techniques Research Articles

Positivity and dynamics preserving discretization schemes for nonlinear evolution equations

Discretization of a continuous-time system of dierential equations becomes inevitable due to the lack of analytical solutions. Standard discretization techniques, however, have many things that could be improved, e.g., the positivity of the solution and dynamic consistency may be lost, and stability and convergence may depend on the step length. A nonstandard nite dierence (NSFD) scheme is sometimes used to avoid inconsistencies. There are two fundamental issues regarding the construction of NSFD models. First, how to construct the denominator function of the discrete rst-order derivative? Second, how to discretize the nonlinear terms of a given dierential equation with nonlocal terms? We dene here a uniform technique for nonlocal discretization and construction of denominator function for NSFD models. We have discretized a couple of highly nonlinear continuous-time population models using these consistent rules. We give analytical proof in each case to show that the proposed NSFD model has identical dynamic properties to the continuous-time model. It is also shown that each NSFD system is positively invariant, and its dynamics do not depend on the step size. Numerical experiments have also been performed in favour of such claims.

SIMPLIFIED ADAPTIVE TRIANGULATION OF THE CONTACT BOUNDARIES OF THE DAM MODEL

To numerically solve the system of integral equations, it is customary to establish a discrete grid within each integration area. In the context of 3D modeling, these areas correspond to surfaces situated in space. The standard discretization technique employed for the computational domain is triangulation. This study addresses the integral equation system pertinent to the electrical tomography of dams. The structural model encompasses an embankment dam, the upstream and downstream water bodies, the dam's base, and a potential leakage region on the upstream side. An alternative configuration may be encountered in specific scenarios with no water downstream. Consequently, the model may incorporate up to nine distinct contact boundaries. Accordingly, the system of integral equations comprises an equivalent number of equations. Effectively resolving this system through numerical methods necessitates applying triangulation techniques to these diverse surfaces. While mathematical packages like Matlab offer triangulation functions, they may not fully address the specific demands of the problem. Additionally, the grid resolution should be heightened in proximity to key elements such as the sounding line, the supply electrode, and the various contact lines within the medium. These considerations transform the triangulation task into a distinct subtask within the numerical simulation of the resistivity tomography problem. In this paper, we provide our specific approach to this problem. The simplification of the triangulation algorithm is rooted in the predominant utilization of the two-dimensional geometric properties inherent to the object under study. For most contact boundaries, the triangulation is constructed layer by layer with a gradual modulation in triangle dimensions as one progresses from one layer to the next, orthogonal to the axis of the dam. Concerning the surface corresponding to the leakage area, cylindrical coordinates are used for surface parameterization. This approach enables partitioning the surface into discrete strata, facilitating a systematic, layer-by-layer grid construction. Additionally, points at the intersections of contact boundaries are integrated into the pre-existing triangulation by applying a standard function within the Matlab package. In the future, the mathematical modeling based on the Integral Equation Method with adaptive discretization will help incorporate real-time computations into information systems related to monitoring hydraulic structures.

Distributed ℋ 2 -Matrices for Boundary Element Methods

Standard discretization techniques for boundary integral equations, e.g., the Galerkin boundary element method, lead to large densely populated matrices that require fast and efficient compression techniques like the fast multipole method or hierarchical matrices. If the underlying mesh is very large, running the corresponding algorithms on a distributed computer is attractive, e.g., since distributed computers frequently are cost-effective and offer a high accumulated memory bandwidth. Compared to the closely related particle methods, for which distributed algorithms are well-established, the Galerkin discretization poses a challenge, since the supports of the basis functions influence the block structure of the matrix and therefore the flow of data in the corresponding algorithms. This article introduces distributed ℋ 2 -matrices, a class of hierarchical matrices that is closely related to fast multipole methods and particularly well-suited for distributed computing. While earlier efforts required the global tree structure of the ℋ 2 -matrix to be stored in every node of the distributed system, the new approach needs only local multilevel information that can be obtained via a simple distributed algorithm, allowing us to scale to significantly larger systems. Experiments show that this approach can handle very large meshes with more than 130 million triangles efficiently.

Modeling curvature-resisting material surfaces with isogeometric analysis

Improved understanding of solid surface energy and its role in the overall mechanical properties is of great interest due to the emerging manufacturing techniques of nanostructures, coatings, and synthetic/biological bilayer–polymer hybrids. Continuum numerical modeling of surface stresses efficiently incorporates a zero-thickness membrane bonded to a bulk, intrinsically accounting for surface tension and surface elasticity. Compressive surface stresses are not possible in a purely membrane formulation, ignoring the surface flexural resistance. The extension of material surfaces to account for flexural resistance, i.e., the Steigmann–Ogden model, requires spatial derivatives of second order, posing significant challenges to standard discretization techniques. Hence, the effect of surface curvature resistance on the overall mechanical behavior of complex geometries remains elusive. Here, we develop a three-dimensional computational formulation of curvature-dependent surface energetics at finite strains using surface-enriched isogeometric analysis. Coupled with a hyperelastic bulk, bending-resistance of material surfaces furnishes a new physical length scale, i.e., the elastobending length. We quantify the effect of elastobending deformations for several numerical examples involving soft materials with thin coatings and liquid-shell surfaces, capturing budding-like behavior observed at cell membranes. Our results demonstrate a stiffer overall mechanical behavior when material surfaces resist bending deformations and illustrate how curvature effects lead to complex budding deformations at non-zero initial curvature states. The proposed methodology provides a robust computational foundation to help improve our understanding and mechanical characterization of soft solids, nanostructures, and biological membranes at small scales.

Driving bifurcating parametrized nonlinear PDEs by optimal control strategies: application to Navier–Stokes equations with model order reduction

This work deals with optimal control problems as a strategy to drive bifurcating solution of nonlinear parametrized partial differential equations towards a desired branch. Indeed, for these governing equations, multiple solution configurations can arise from the same parametric instance. We thus aim at describing how optimal control allows to change the solution profile and the stability of state solution branches. First of all, a general framework for nonlinear optimal control problem is presented in order to reconstruct each branch of optimal solutions, discussing in detail the stability properties of the obtained controlled solutions. Then, we apply the proposed framework to several optimal control problems governed by bifurcating Navier–Stokes equations in a sudden-expansion channel, describing the qualitative and quantitative effect of the control over a pitchfork bifurcation, and commenting in detail the stability eigenvalue analysis of the controlled state. Finally, we propose reduced order modeling as a tool to efficiently and reliably solve parametric stability analysis of such optimal control systems, which can be challenging to perform with standard discretization techniques such as Finite Element Method.

A separated representation involving multiple time scales within the Proper Generalized Decomposition framework

Solutions of partial differential equations can exhibit multiple time scales. Standard discretization techniques are constrained to capture the finest scale to accurately predict the response of the system. In this paper, we provide an alternative route to circumvent prohibitive meshes arising from the necessity of capturing fine-scale behaviors. The proposed methodology is based on a time-separated representation within the standard Proper Generalized Decomposition, where the time coordinate is transformed into a multi-dimensional time through new separated coordinates, each representing one scale, while continuity is ensured in the scale coupling. For instance, when considering two different time scales, the governing Partial Differential Equation is commuted into a nonlinear system that iterates between the so-called microtime and macrotime, so that the time coordinate can be viewed as a 2D time. The macroscale effects are taken into account by means of a finite element-based macro-discretization, whereas the microscale effects are handled with unidimensional parent spaces that are replicated throughout the time domain. The resulting separated representation allows us a very fine time discretization without impacting the computational efficiency. The proposed formulation is explored and numerically verified on thermal and elastodynamic problems.

Fast Computation of Multi-Parametric Electromagnetic Fields in Synchronous Machines by Using PGD-Based Fully Separated Representations

A novel Model Order Reduction (MOR) technique is developed to compute high-dimensional parametric solutions for electromagnetic fields in synchronous machines. Specifically, the intrusive version of the Proper Generalized Decomposition (PGD) is employed to simulate a Permanent-Magnet Synchronous Motor (PMSM). The result is a virtual chart allowing real-time evaluation of the magnetic vector potential as a function of the operation point of the motor, or even as a function of constructive parameters, such as the remanent flux in permanent magnets. Currently, these solutions are highly demanded by the industry, especially with the recent developments in the Electric Vehicle (EV). In this framework, standard discretization techniques require highly time-consuming simulations when analyzing, for instance, the noise and vibration in electric motors. The proposed approach is able to construct a virtual chart within a few minutes of off-line simulation, thanks to the use of a fully separated representation in which the solution is written from a series of functions of the space and parameters coordinates, with full space separation made possible by the use of an adapted geometrical mapping. Finally, excellent performances are reported when comparing the reduced-order model with the more standard and computationally costly Finite Element solutions.

Impact of Catchment Discretization and Imputed Radiation on Model Response: A Case Study from Central Himalayan Catchment

Distributed and semi-distributed hydrological modeling approaches commonly involve the discretization of a catchment into several modeling elements. Although some modeling studies were conducted using triangulated irregular networks (TINs) previously, little attention has been given to assess the impact of TINs as compared to the standard catchment discretization techniques. Here, we examine how different catchment discretization approaches and radiation forcings influence hydrological simulation results. Three catchment discretization methods, i.e., elevation zones (Hypsograph) (HYP), regular square grid (SqGrid), and TIN, were evaluated in a highly steep and glacierized Marsyangdi-2 river catchment, central Himalaya, Nepal. To evaluate the impact of radiation on model response, shortwave radiation was converted using two approaches: one with the measured solar radiation assuming a horizontal surface and another with a translation to slopes. The results indicate that the catchment discretization has a great impact on simulation results. Evaluation of the simulated streamflow value using Nash–Sutcliffe efficiency (NSE) and log-transformed Nash–Sutcliffe efficiency (LnNSE) shows that highest model performance was obtained when using TIN followed by HYP (during the high flow condition) and SqGrid (during the low flow condition). Similar order of precedence in relative model performance was obtained both during the calibration and validation periods. Snow simulated from the TIN-based discretized models was validated with Moderate Resolution Imaging Spectroradiometer (MODIS) snow products. Critical Success Indexes (CSI) between TIN-based discretized model snow simulation and MODIS snow were found satisfactory. Bias in catchment average snow cover area from the models with and without using imputed radiation is less than two percent, but implementation of imputed radiation into the Statkraft Hydrological Forecasting Toolbox (Shyft) gives better CSI with MODIS snow.

Multiscale proper generalized decomposition based on the partition of unity

SummarySolutions of partial differential equations could exhibit a multiscale behavior. Standard discretization techniques are constraints to mesh up to the finest scale to predict accurately the response of the system. The proposed methodology is based on the standard proper generalized decomposition rationale; thus, the PDE is transformed into a nonlinear system that iterates between microscale and macroscale states, where the time coordinate could be viewed as a 2D time, representing the microtime and macrotime scales. The macroscale effects are taken into account because of an FEM‐based macrodiscretization, whereas the microscale effects are handled with unidimensional parent spaces that are replicated throughout the domain. The proposed methodology can be seen as an alternative route to circumvent prohibitive meshes arising from the necessity of capturing fine‐scale behaviors.

On the dynamic consistency of a two-species competitive discrete system with toxicity: Local and global analysis

Discretization of a highly nonlinear system is inevitable for good approximation of its solutions due to lack of its analytical solution. However, a discrete model constructed by standard discretization technique frequently produces spurious dynamics which are not observed in its constituent continuous model. It is therefore important to construct a discrete model that produces the exact dynamics of its continuous counterpart. Here we transform a continuous-time two species nonlinear competition model with the effect of toxic substances into a discrete-time system to show its dynamic consistencies with the continuous system by employing nonstandard finite difference scheme. Linear stability property of each fixed point of the derived model has been proven to be identical with those of the continuous system and the dynamics has been shown to be independent of the step-size. Moreover, we prove the global stability of the coexisting equilibrium point of our discrete model. Comparative simulation results demonstrate that the standard discrete system shows complex dynamics, like flip bifurcation and chaos depending on the step-size, but the proposed discrete model does not show any such spurious dynamics.

Empirical Interpolation Decomposition

Many physical problems need a multidimensional description and involve high dimensional spaces. Standard discretization techniques often lead to an excessive computation time. To solve this problem, we develop in this paper an empirical interpolation decomposition (EID) for multivariate functions. This method provides an approximate representation of a given function in separate form. Error estimates of the developed EID are derived and some properties are given.

English

The paper focuses on a low-rank tensor structured representation of Slatertype and Hydrogen-like orbital basis functions that can be used in electronic structure calculations. Standard packages use the Gaussian-type basis functions which allow us to analytically evaluate the necessary integrals. Slater-type and Hydrogen-like orbital functions are physically more appropriate, but they are not analytically integrable. A numerical integration is too expensive when using the standard discretization techniques due the dimensionality of the problem. However, it can be effectively performed using the tensor representation of basis functions. Furthermore, this approach can take advantage of parallel computing.

Preconditioned Low-rank Riemannian Optimization for Linear Systems with Tensor Product Structure

The numerical solution of partial differential equations on high-dimensional domains gives rise to computationally challenging linear systems. When using standard discretization techniques, the size of the linear system grows exponentially with the number of dimensions, making the use of classic iterative solvers infeasible. During the last few years, low-rank tensor approaches have been developed that allow one to mitigate this curse of dimensionality by exploiting the underlying structure of the linear operator. In this work, we focus on tensors represented in the Tucker and tensor train formats. We propose two preconditioned gradient methods on the corresponding low-rank tensor manifolds: a Riemannian version of the preconditioned Richardson method as well as an approximate Newton scheme based on the Riemannian Hessian. For the latter, considerable attention is given to the efficient solution of the resulting Newton equation. In numerical experiments, we compare the efficiency of our Riemannian algorit...

Efficient Stabilization of Advection Terms Involved in Separated Representations of Boltzmann and Fokker-Planck Equations

AbstractThe fine description of complex fluids can be carried out by describing the evolution of each individual constituent (e.g. each particle, each macromolecule, etc.). This procedure, despite its conceptual simplicity, involves many numerical issues, the most challenging one being that related to the computing time required to update the system configuration by describing all the interactions between the different individuals. Coarse grained approaches allow alleviating the just referred issue: the system is described by a distribution function providing the fraction of entities that at certain time and position have a particular conformation. Thus, mesoscale models involve many different coordinates, standard space and time, and different conformational coordinates whose number and nature depend on the particular system considered. Balance equation describing the evolution of such distribution function consists of an advection-diffusion partial differential equation defined in a high dimensional space. Standard mesh-based discretization techniques fail at solving high-dimensional models because of the curse of dimensionality. Recently the authors proposed an alternative route based on the use of separated representations. However, until now these approaches were unable to address the case of advection dominated models due to stabilization issues. In this paper this issue is revisited and efficient procedures for stabilizing the advection operators involved in the Boltzmann and Fokker-Planck equation within the PGD framework are proposed.

Separated representations of 3D elastic solutions in shell geometries

Abstract Background The solution of 3D models in degenerated geometries in which some characteristic dimensions are much lower than the other ones -e.g. beams, plates, shells,...- is a tricky issue when using standard mesh-based discretization techniques. Methods Separated representations allow decoupling the meshes used for approximating the solution along each coordinate. Thus, in plate or shell geometries 3D solutions can be obtained from a sequence of 2D and 1D problems allowing fine and accurate representation of the solution evolution along the thickness coordinate while keeping the computational complexity characteristic of 2D simulations. In a former work this technique was considered for addressing the 3D solution of thermoelastic problems defined in plate geometries. In this work, the technique is extended for addressing the solution of 3D elastic problems defined in shell geometries. Results The capabilities of the proposed approach are illustrated by considering some numerical examples involving different degrees of complexity, from simple shells to composite laminates involving stiffeners. Conclusions The analyzed examples prove the potentiality and efficiency of the proposed strategy, where the computational complexity was found evolving as reported in our former works, proving that 3D solutions can be computed at a 2D cost.

Calculation of the interface curvature and normal vector with the level-set method

This article addresses the use of the level-set method for capturing the interface between two fluids. One of the advantages of the level-set method is that the curvature and the normal vector of the interface can be readily calculated from the level-set function. However, in cases where the level-set method is used to capture topological changes, the standard discretization techniques for the curvature and the normal vector do not work properly. This is because they are affected by the discontinuities of the signed-distance function half-way between two interfaces. This article addresses the calculation of normal vectors and curvatures with the level-set method for such cases. It presents a discretization scheme based on the geometry-aware curvature discretization by Macklin and Lowengrub [1]. As the present scheme is independent of the ghost-fluid method, it becomes more generally applicable, and it can be implemented into an existing level-set code more easily than Macklin and Lowengrub’s scheme [1]. The present scheme is compared with the second-order central-difference scheme and with Macklin and Lowengrub’s scheme [1], first for a case with no flow, then for a case where two drops collide in a 2D shear flow, and finally for a case where two drops collide in an axisymmetric flow. In the latter two cases, the Navier–Stokes equations for incompressible two-phase flow are solved. The article also gives a comparison of the calculation of normal vectors with the direction difference scheme presented by Macklin and Lowengrub in [2] and with the present discretization scheme. The results show that the present discretization scheme yields more robust calculations of the curvature than the second-order central difference scheme in areas where topological changes are imminent. The present scheme compares well to Macklin and Lowengrub’s method [1]. The results also demonstrate that the direction difference scheme [2] is not always sufficient to accurately calculate the normal vectors.

Robust Solution of Singularly Perturbed Problems Using Multigrid Methods

We consider the problem of solving linear systems of equations that arise in the numerical solution of singularly perturbed ordinary and partial differential equations of reaction-diffusion type. Standard discretization techniques are not suitable for such problems and, so, specially tailored methods are required, usually involving adapted or fitted meshes that resolve important features such as boundary and/or interior layers. In this study, we consider classical finite difference schemes on the layer adapted meshes of Shishkin and Bakhvalov. We show that standard direct solvers exhibit poor scaling behavior, with respect to the perturbation parameter, when solving the resulting linear systems. We propose and prove optimality of a new block-structured preconditioning approach that is robust for small values of the perturbation parameter, and compares favorably with standard robust multigrid preconditioners for these linear systems. We also derive stopping criteria which ensure that the potential accuracy of the layer-resolving meshes is achieved.

Effective models for CO2 migration in geological systems with varying topography

Geological CO2 sequestration relies on a competent sealing layer, or caprock, that bounds the formation top and prevents vertical migration of CO2 and brine. Modeling studies have shown that caprock topography, or roughness, can significantly decrease updip migration speed of CO2 and increase structural trapping. Caprock roughness can be characterized at different spatial scales. For instance, regional‐scale features such as domes, traps, and spill points can be detected in seismic surveys and have been shown to affect large‐scale migration patterns. However, structural and topographical variability, known as rugosity, exists below seismic detection limits but can be measured at the scale of centimeters and meters using LiDAR scanning of formation outcrops. Little is known about the actual impact of structural rugosity on CO2 plume migration. Practically speaking, given the large scales required to model commercial scale CO2 storage projects and the limitations on computational power, only seismic‐scale caprock topography can be resolved using standard discretization techniques. Therefore, caprock variability that exists below the model resolution scale is defined as subscale and must be handled by upscaling. In this paper we derive effective equations for CO2 migration that include the impact of fine‐scale variability in caprock topography using static equilibrium upscaling, an approach that is adapted for the vertical equilibrium modeling framework. The effective equations give estimates of the impact of rugosity on CO2 plume migration and trapping in large‐scale systems.

A fully discrete Galerkin method for high frequency exterior acoustic scattering in three dimensions

Standard Galerkin discretization techniques (with locally- or globally-supported basis functions) for boundary integral equations are inefficient for high frequency three dimensional exterior scattering simulations because they require a fixed number of unknowns per wavelength in each dimension, leading to large CPU time and memory requirements to set up the dense Galerkin matrix, with each entry requiring evaluation of multi-dimensional highly oscillatory integrals. In this work, using globally-supported basis functions, we describe an efficient fully discrete Galerkin surface integral equation algorithm for simulating high frequency acoustic scattering by three dimensional convex obstacles that includes a powerful integration scheme for evaluation of four dimensional Galerkin integrals with high-order accuracy. Such high-order order accuracy for various practically relevant frequencies ( k ∈ [1, 100,000]) substantially improves on approximations based on standard asymptotic techniques. We demonstrate the efficiency of our algorithm for spherical and non-spherical convex scattering for several wavenumbers 1 ⩽ k ⩽ 100,000 for low to high order prescribed tolerance. Our fully discrete algorithm requires only mild growth in the number of unknowns and CPU time as the frequency increases.

Efficient mold cooling optimization by using model reduction

Optimization and inverse identification are two procedures usually encountered in many industrial processes reputed gourmand for the computing time view point. In fact, optimization implies to propose a trial solution whose accuracy is then evaluated, and if needed it must be updated in order to minimize a certain cost function. In the case of mold cooling optimization the evaluation of the solution quality needs the solution of a thermal model, in the whole domain and during the thermal history. Thus, the optimization process needs several iterations and then the computational cost can become enormous. In this work we propose the use of model reduction for accomplishing this kind of simulations. Thus, only one thermal model is solved using the standard discretization technique. After that, the most important modes defining the temperature evolution are extracted by invoking the proper orthogonal decomposition, and all the other thermal model solutions are performed by using the reduced order approximation basis just extracted. The CPU time savings can be impressive.