Integral Solver Research Articles

The impact of non-monotonic relaxation on numerical accuracy and stability of viscoelastic finite elements

Analysis of sound and vibration in linear viscoelastic materials is based on a relaxation function that is the time-dependent stress resulting from a unit step in strain. It is well established that relaxation functions decrease monotonically in time. The most common constitutive relaxation functions, such as the Kelvin–Voigt or General Maxwell Model, produce curves that must be monotonic. However, recently, the authors have produced a methodology that allows the use of discrete data for the time-dependent relaxation function. One result of this approach is the possibility of non-monotonic relaxation function input. This non-monotonic data may be the result of test and measurement errors. Alternatively, this non-monotonicity may be an actual anomalous physical result, such as in rock salt (He et al. 2019). The authors use validated viscoelastic finite element models in the time domain to determine the instability caused by such non-monotonic datasets. The study uses randomized non-monotonic datasets that maintain the fundamental material trend to find unstable or inaccurate results. A mathematical review highlights the conditions causing numerical instability in modern stepwise time integration solvers. The results provide clarity for viscoelastic analysts considering low-precision viscoelastic measurements or unusual physical properties. [Work supported by ONR under Grant N00014-22-1-2785.]

Optimized parallelization of boundary integral Poisson-Boltzmann solvers

The Poisson-Boltzmann (PB) model governs the electrostatics of solvated biomolecules, i.e., potential, field, energy, and force. These quantities can provide useful information about protein properties, functions, and dynamics. By considering the advantages of current algorithms and computer hardware, we focus on the parallelization of the treecode-accelerated boundary integral (TABI) PB solver using the Message Passing Interface (MPI) on CPUs and the direct-sum boundary integral (DSBI) PB solver using KOKKOS on GPUs. We provide optimization guidance for users when the DSBI solver on GPU or the TABI solver with MPI on CPU should be used depending on the size of the problem. Specifically, when the number of unknowns is smaller than a predetermined threshold, the GPU-accelerated DSBI solver converges rapidly thus has the potential to perform PB model-based molecular dynamics or Monte Carlo simulation. As practical applications, our parallelized boundary integral PB solvers are used to solve electrostatics on selected proteins that play significant roles in the spread, treatment, and prevention of COVID-19 virus diseases. For each selected protein, the simulation produces the electrostatic solvation energy as a global measurement and electrostatic surface potential for local details.

Faster diffusion model with improved quality for particle cloud generation

Building on the success of PC-JeDi we introduce PC-Droid, a substantially improved diffusion model for the generation of jet particle clouds. By leveraging a new diffusion formulation, studying more recent integration solvers, and training on all jet types simultaneously, we are able to achieve state-of-the-art performance for all types of jets across all evaluation metrics. We study the trade-off between generation speed and quality by comparing two attention based architectures, as well as the potential of consistency distillation to reduce the number of diffusion steps. Both the faster architecture and consistency models demonstrate performance surpassing many competing models, with generation time up to two orders of magnitude faster than PC-JeDi and three orders of magnitude faster than . Published by the American Physical Society 2024

Receding Galerkin Optimal Control with High-Order Sliding Mode Disturbance Observer for a Boiler-Turbine Unit

The control of the boiler-turbine unit is important for its sustainable and robust operation in power plants, which faces great challenges due to the control unit’s serious nonlinearity, unmeasurable states, variable constraints, and unknown time-varying lumped disturbances. To address the above issues, this paper proposes a receding Galerkin optimal controller with a high-order sliding mode disturbance observer in a composite scheme, in which a high-order sliding mode disturbance observer is first employed to estimate the lumped disturbances based on a deviation form of the mathematical model of the boiler-turbine unit. Subsequently, under the hypothesis of state constraint, a receding Galerkin optimal controller is designed to compensate the lumped disturbances by embedding their estimates into the mathematically based predictive model at each sampling time instant. With the help of an interpolation polynomial, Gauss integration, and nonlinear solvers, an optimal control law is then obtained based on a Galerkin optimization algorithm. Consequently, disturbance rejection, target tracking, and constraint handling performance of a controlled closed-loop system are improved. Some simulation cases are conducted on a mathematical boiler-turbine unit model to demonstrate the effectiveness of the proposed method, which is supported by the quantitative result analysis, such as tracking and disturbance rejection performance indexes.

ARKODE: A Flexible IVP Solver Infrastructure for One-step Methods

We describe the ARKODE library of one-step time integration methods for ordinary differential equation (ODE) initial-value problems (IVPs). In addition to providing standard explicit and diagonally implicit Runge–Kutta methods, ARKODE supports one-step methods designed to treat additive splittings of the IVP, including implicit-explicit (ImEx) additive Runge–Kutta methods and multirate infinitesimal (MRI) methods. We present the role of ARKODE within the SUNDIALS suite of time integration and nonlinear solver libraries, the core ARKODE infrastructure for utilities common to large classes of one-step methods, as well as its use of “time stepper” modules enabling easy incorporation of novel algorithms into the library. Numerical results show example problems of increasing complexity, highlighting the algorithmic flexibility afforded through this infrastructure, and include a larger multiphysics application leveraging multiple algorithmic features from ARKODE and SUNDIALS.

An integral equation method for the advection-diffusion equation on time-dependent domains in the plane

Boundary integral methods are attractive for solving homogeneous elliptic partial differential equations on complex geometries, since they can offer accurate solutions with a computational cost that is linear or close to linear in the number of discretization points on the boundary of the domain. However, these numerical methods are not straightforward to apply to time-dependent equations, which often arise in science and engineering. We address this problem with an integral equation-based solver for the advection-diffusion equation on moving and deforming geometries in two space dimensions. In this method, an adaptive high-order accurate time-stepping scheme based on semi-implicit spectral deferred correction is applied. One time-step then involves solving a sequence of non-homogeneous modified Helmholtz equations, a method known as elliptic marching. Our solution methodology utilizes several recently developed methods, including special purpose quadrature, a function extension technique and a spectral Ewald method for the modified Helmholtz kernel. Special care is also taken to handle the time-dependent geometries. The numerical method is tested through several numerical examples to demonstrate robustness, flexibility and accuracy.

An accelerated subspaces recycling strategy for the deflation of parametric linear systems based on model order reduction

Krylov subspace recycling has been extensively used to facilitate the solution of sequences of linear systems by constructing a deflation subspace and accelerating the convergence of a corresponding iterative solver. However, most existing techniques update the recycled subspace sequentially for each system, thus inducing a potentially high computational cost. In that context, this work proposes a method to accelerate the above procedure for the case of multiresolution analyses of affine parametric systems, by decoupling the solution of the system from the construction of the recycled subspace. The proposed method follows the projection based Model Order Reduction (MOR) logic that splits operations into an offline and an online phase and therefore proves particularly beneficial in case of affine parametrizations that involve multiple affine coefficients. In the offline phase of the method an especially tailored version of the Automatic Krylov subspaces Recycling (AKR) algorithm (Panagiotopoulos et al., 2021), proposed within this work and denoted as AKR-D, is employed. In brief, AKR-D constructs a high quality recycling basis W for a predefined parameter interval Ψ by targeting a desirable convergence rate at an automatically generated set of parameter values Ω⊂Ψ. The upfront construction of W enables an a-priori Galerkin projection of the affinely described full parametric system to yield a reduced order model (ROM). This ROM can be leveraged in the online phase of the proposed recycling method to accelerate the construction of the corresponding projector P onto W=span{W}, whose implicit assembly is required in most recycling based deflation schemes. The effectiveness of the proposed strategy is investigated in terms of algorithmic complexity and is demonstrated on a densely parametrized system arising within the boundary integral solver of the Helmholtz equation and a random sparse parametrization example.

TABI-PB 2.0: An Improved Version of the Treecode-Accelerated Boundary Integral Poisson-Boltzmann Solver.

This work describes TABI-PB 2.0, an improved version of the treecode-accelerated boundary integral Poisson-Boltzmann solver. The code computes the electrostatic potential on the molecular surface of a solvated biomolecule, and further processing yields the electrostatic solvation energy. The new implementation utilizes the NanoShaper surface triangulation code, node-patch boundary integral discretization, a block preconditioner, and a fast multipole method based on barycentric Lagrange interpolation and dual tree traversal. Performance-critical portions of the code were implemented on a GPU. Numerical results for protein 1A63 and two viral capsids (Zika, H1N1) demonstrate the code's accuracy and efficiency.

Finite Element Modeling of Atmospheric Water Extraction by Way of Highly Porous Adsorbents: A Roadmap for Solver Construction with Model Factor Sensitivity Screening.

A finite element model (FEM) is developed for use in determining adsorption system performance. The model is intended to guide novel adsorbent structure fabrication and atmospheric water harvesting device design. We survey a variety of governing equation factor inputs and relationships which describe the interaction between zeolite 13X and water vapor. Mitigation strategies are discussed for detecting the breakdown of continuum modeling at the microscale wherein Knudsen effects and other anomalous behaviors emerge. Characterization of model factor inputs and the techniques for their sourcing is described with consideration to the construction of a high throughput multiscale shape optimized computational schema. Four objectives guided the development of this model. Our first objective was to understand the implementation of adsorption system equations and the assumptions that could prevent reliable predictability. The second objective was to assemble, reduce, and analyze model constants and approximations that express FEM coefficient calculations as physical forces and thermodynamic properties which could be derived from other computational methods. Third, we analyzed factor sensitivity of model inputs by way of a 2k factorial screening to determine which inputs are driving the physics of water harvesting adsorption systems. The fourth objective was to design the FEM solver for integration into a multiscale high throughput topologically optimized schema. The main finding of the solver factor screening indicates that total micropore volume has the highest value characteristics in relation to water uptake.

Approximate weighted model integration on DNF structures

Weighted model counting consists of computing the weighted sum of all satisfying assignments of a propositional formula. Weighted model counting is well-known to be #P-hard for exact solving, but admits a fully polynomial randomized approximation scheme when restricted to DNF structures. In this work, we study weighted model integration, a generalization of weighted model counting which involves real variables in addition to propositional variables, and pose the following question: Does weighted model integration on DNF structures admit a fully polynomial randomized approximation scheme? Building on classical results from approximate weighted model counting and approximate volume computation, we show that weighted model integration on DNF structures can indeed be approximated for a class of weight functions. Our approximation algorithm is based on three subroutines, each of which can be a weak (i.e., approximate), or a strong (i.e., exact) oracle, and in all cases, comes along with accuracy guarantees. We experimentally verify our approach over randomly generated DNF instances of varying sizes, and show that our algorithm scales to large problem instances, involving up to 1K variables, which are currently out of reach for existing, general-purpose weighted model integration solvers.

Power performance comparison analysis of FPGA using Riemann Definite Integral Equation Solver

Definite integral equation solvers were hitherto developed using software which take long time to generate results. The complexity of designs and verifications, due to advancement in calculus and linear algebra, data dependency in designs, and the need for parallelism in models, necessitated the need for faster, accurate, and reliable hardware-based solvers. This will reduce latency and excessive memory requirement in generating values of definite integrals. The Field Programmable Gate Array (FPGA)-Based Definite Integral Solver developed in this research is able to solve definite integral equations with less resource utilization (1.08% look up tables (LUT), 0.07% flip flops (FF), 10.81% digital signal processing (DSP) and 23.86% input-output (IO)) and hence minimum switching activities at lower power (Signals: 0.421W, Logic: 0.345W, Digital Signal Processing: 1.904W and Input/Output: 0.015W) that translates to fast and accurate results generation. The power and resource utilization evaluation on the solver and the performance comparison with similar works showed that excessive memory consumption and resource utilization has been reduced significantly by the design.

Foundry-fabricated grating coupler demultiplexer inverse-designed via fast integral methods

Silicon photonics is an emerging technology which, enabling nanoscale manipulation of light on chips, impacts areas as diverse as communications, computing, and sensing. Wavelength division multiplexing is commonly used to maximize throughput over a single optical channel by modulating multiple data streams on different wavelengths concurrently. Traditionally, wavelength (de)multiplexers are implemented as monolithic devices, separate from the grating coupler, used to couple light into the chip. This paper describes the design and measurement of a grating coupler demultiplexer—a single device which combines both light coupling and demultiplexing capabilities. The device was designed by means of a custom inverse design algorithm which leverages boundary integral Maxwell solvers of extremely rapid convergence as the mesh is refined. To the best of our knowledge, the fabricated device enjoys the lowest insertion loss reported for grating demultiplexers, small size, high splitting ratio, and low coupling-efficiency imbalance between ports, while meeting the fabricability constraints of a standard UV lithography process.

Nested kernel degeneration-based boundary element method solver for rapid computation of eddy current signals

In this article, the nested kernel degeneration (NKD) based integral equation fast solver is proposed to compress the dense system matrix generated by boundary element method (BEM) for calculating eddy current signals in 3D nondestructive evaluation (NDE) problems. The near and diagonal block interactions are computed and stored by full matrices. However, for the far block interactions, not like in the single level kernel degeneration-based method that all the interaction matrices are compressed by the interpolations, the nested approximation framework makes the expression of the higher-level matrices by the ones at leaf levels and the transfer matrices between neighboring levels. Different integral kernels are degenerated and nested to ensure the performance of NKD based method. The robustness and efficiency of the proposed solver are validated and verified by making numerical predictions and comparing them with the ones achieved by analytical method, numerical methods and experiments. Our numerical experiments show that the NKD solver is more efficient than the butterfly solver (multilevel adaptive cross approximation algorithm (MLACA)) for electrically small eddy current NDE problem. With the NKD based method, O(N) computational complexity is achieved for multiscale low frequency eddy current simulations without sacrificing accuracy, where N is the number of unknowns. So far, to our best knowledge, the proposed method with linear complexity is the most efficient integral equation-based BEM forward solver to predict the eddy current signals from cracks and flaws.

A Coupled Inviscid–Viscous Airfoil Analysis Solver, Revisited

This paper presents a self-contained, comprehensive exposition and new elements of a classic airfoil analysis technique: an integral boundary-layer solver coupled to a vortex-panel potential-flow method. The resulting solver MFOIL, implemented in MATLAB and made freely available, builds on the XFOIL code and documentation, which serve as the inspiration, reference, and standard of comparison. Modifications made in the present implementation include an augmented-state coupled solver for more control in limiting the state update, a new stagnation-point formulation to reduce leading-edge oscillations in the boundary-layer variables, and a more robust treatment of the amplification rate near transition. Several results highlight these modifications and demonstrate capabilities of the code.

SMIS: A Stepwise Multiple Integration Solver Using a CAS

Multiple Integration is a very important topic in different applications in Engineering and other Sciences. Using numerical software to get an approximation to the solution is a normal procedure. Another approach is working in an algebraic form to obtain an exact solution or to get general solutions depending on different parameters. Computer Algebra Systems (CAS) are needed for this last approach. In this paper, we introduce SMIS, a new stepwise solver for multiple integration developed in a CAS. The two main objectives of SMIS are: (1) to increase the capabilities of CAS to help the user to deal with this topic and (2) to be used in Math Education providing an important tool for helping with the teaching and learning process of this topic. SMIS can provide just the final solution or an optional stepwise solution (even including some theoretical comments). The optional stepwise solutions provided by SMIS are of great help for (2). Although SMIS has been developed in the specific CAS Derive, since the code is provided, it can be easily migrated to any CAS which deals with integrals and text management that allow us to display comments for intermediate steps.

Material and geometric effects on propulsion of a fish tail

We investigate the effects of material flexibility and aspect ratio on the propulsion of flapping tails. The tail, which is assumed to deform in the bending direction only, is modeled using the Euler–Bernoulli beam theory. The hydrodynamic loads generated by the flapping motion are calculated using the three-dimensional unsteady vortex lattice method. The finite element method is used to solve the coupled time-dependent equations of motion using an implicit solver for time integration. The results show improvement in the thrust and propulsive efficiency over a specific range of non-dimensional flexibility defined by the ratio of the elastic forces to fluid pressure forces. Structural and flow characteristics associated with the improved performance are discussed. As for geometric effects, the performance depends on the excitation frequency. At low frequencies, the improvement is continuous with increasing the aspect ratio in a manner similar to that of rigid tails. At higher frequencies, the improvement is limited to a region defined by aspect ratios that are less than 0.5. The extent of the improvement depends on the flexibility.

Integral equation fast solver with truncated and degenerated kernel for computing flaw signals in eddy current non-destructive testing

In this article, the kernel independent H-matrix method is applied as the mathematical framework. The dense matrix generated by the selected integral equation, without low frequency breakdown problem for solving 3D arbitrary shaped eddy current nondestructive evaluation (NDE) forward problems, is compressed by the truncated and degenerated (TD) kernel function method. It results in a significant reduction in computational burden with nearly linear complexity. The kernel functions (Green's function) are degenerated by Lagrange polynomials which leads to two advantages: firstly, the double integrals are separated in two single integrals, and secondly, they are represented in a factorized form with good accuracy control. Meanwhile, due to the nature of the kernel function with exponential decay in the lossy medium, it is truncated to improve the overall performance. The TD kernel function-based boundary element method (BEM) fits well for solving eddy current NDE problems over the adaptive cross approximation based BEM, especially for detections of flaws or slots in planar objects. Several numerical predictions calculated by the proposed method are compared with those achieved by others such as the experiment, analytical and semi-analytical methods from benchmarks. The robustness and efficiency are demonstrated with excellent agreements and reduced complexity.

Relationship between the contact force strength and numerical inaccuracies in piecewise-smooth systems

This work studies the different types of behavior and inaccuracies that can occur when contact is not adequately accounted for in a dynamical system with freeplay, as the strength of the contact stiffness increases. The MATLAB® ode45 time integration solver, with the built-in Event Location capability, is first validated using past experimental data from a forced Duffing oscillator with freeplay. Next, numerical results utilizing event location are compared to results neglecting event location in order to highlight possible numerical errors and effects on multistable dynamical responses. Inaccuracies tend to occur in two different ways. First, neglecting event location can affect the boundaries between basins of attraction. Second, neglecting event location has little effect on the behaviors of the attractor solutions themselves besides merely resembling poorly converged solutions. Errors are less pronounced at the limits of soft or hard contact stiffness. This study shows the importance of accurately solving piecewise-smooth systems and the existing correlation between the strength of the contact force and possible numerical inaccuracies.

Fully implicit spectral boundary integral computation of red blood cell flow

This paper is on an implicit time integration scheme for simulation of red blood cell (RBC) flow in an ambient fluid. The intra- and extracellular plasmas are modeled as Stokes flows and represented by boundary integral equations (BIE) written in a weakly singular form. The cell membrane is modeled as a thin elastic shell. Expressed in this way, the RBC flow model constitutes an implicit ordinary differential equation (IODE) in the cell shape. The cell shape and velocity field are discretized spatially by a spectral approach using spherical harmonic basis functions. It is then convenient to express the BIE in the Galerkin form with the spherical harmonics themselves as test functions. The key aspect in this paper is the recognition of the IODE structure of the RBC flow model and consequent application of a multi-step implicit solver for time integration. As with any implicit solver, a nonlinear equation in the cell shape is solved at each time step, for which Newton's method is applied. This requires the Jacobian of the IODE, or equivalently computation of Jacobian-vector products. An important contribution is the formulation of such Jacobian-vector products as evaluating a second BIE. The original weakly singular form is crucial in facilitating this formulation. The implicit solver employs variable order and adaptive time stepping controlled by truncation error and convergence of Newton iterations. Numerical examples show that larger time steps are possible and that the number of matrix-vector products is comparable to explicit methods. Source code is provided in the online supplementary material.

On the breaking inception of unsteady water wave packets evolving in the presence of constant vorticity

Abstract