Operator Splitting Technique Research Articles

Decoupled and energy stable schemes for phase-field surfactant model based on mobility operator splitting technique

In this paper, we investigate numerical methods for the phase-field surfactant (PFS) model, which is a gradient flow system consisting of two nonlinearly coupled Cahn-Hilliard type equations. The main challenge in developing high-order efficient energy stable methods for this system results from the nonlinearity and the strong coupling in the two variables in the free energy functional. We propose two fully decoupled, linear and energy stable schemes based on a linear stabilization approach and an operator splitting technique. We rigorously prove that both schemes can preserve the original energy dissipation law. The techniques employed in these schemes are then summarized into an innovative approach, which we call the mobility operator splitting (MOS), to design high-order decoupled energy stable schemes for a wide class of gradient flow systems. As a particular case, MOS allows different time steps for updating respective variables, leading to a multiple time-stepping strategy for fast-slow dynamics and thus serious improvement of computational efficiency. Various numerical experiments are presented to validate the accuracy, efficiency and other desired properties of the proposed schemes. In particular, detailed phenomena in thin-film pinch-off dynamics can be clearly captured by using the proposed schemes.

Efficient, multimodal, and derivative-free bayesian inference with Fisher–Rao gradient flows

Abstract In this paper, we study efficient approximate sampling for probability distributions known up to normalization constants. We specifically focus on a problem class arising in Bayesian inference for large-scale inverse problems in science and engineering applications. The computational challenges we address with the proposed methodology are: (i) the need for repeated evaluations of expensive forward models; (ii) the potential existence of multiple modes; and (iii) the fact that gradient of, or adjoint solver for, the forward model might not be feasible.

While existing Bayesian inference methods meet some of these challenges individually, we propose a framework that tackles all three systematically. Our approach builds upon the Fisher-Rao gradient flow in probability space, yielding a dynamical system for probability densities that converges towards the target distribution at a uniform exponential rate. This rapid convergence is advantageous for the computational burden outlined in (i). We apply Gaussian mixture approximations with operator splitting techniques to simulate the flow numerically; the resulting approximation can capture multiple modes thus addressing (ii). Furthermore, we employ the Kalman methodology to facilitate a derivative-free update of these Gaussian components and their respective weights, addressing the issue in (iii).

The proposed methodology results in an efficient derivative-free posterior approximation method, flexible enough to handle multi-modal distributions: Gaussian Mixture Kalman Inversion (GMKI). The effectiveness of GMKI is demonstrated both theoretically and numerically in several experiments with multimodal target distributions, including proof-of-concept and two-dimensional examples, as well as a large-scale application: recovering the Navier-Stokes initial condition from solution data at positive times.



Full non-LTE spectral line formation

In the present paper we consider the full nonlocal thermodynamic equilibrium (non-LTE) radiation transfer problem. This formalism allows us to account for deviation from equilibrium distribution of both the radiation field and the massive particles. In the present study, two-level atoms with broadened upper level represent the massive particles. In the absence of velocity-changing collisions, we demonstrate the analytic equivalence of the full non-LTE source function with the corresponding standard non-LTE partial frequency redistribution (PFR) model. We present an iterative method based on operator splitting techniques that can be used to numerically solve the problem at hand. We benchmark it against the standard non-LTE transfer problem for a two-level atom with PFR. We illustrate the deviation of the velocity distribution function of excited atoms from the equilibrium distribution. We also discuss the dependence of the emission profile and the velocity distribution function on elastic collisions and velocity-changing collisions.

PhoenixMR: A GPU-based MRI simulation framework with runtime-dynamic code execution.

Simulations of physical processes and behavior can provide unique insights and understanding of real-world problems. Magnetic Resonance Imaging (MRI) is an imaging technique with several components of complexity. Several of these components have been characterized and simulated in the past. However, several computational challenges prevent simulations from being simultaneously fast, flexible, and accurate. The simulation of MRI experiments is underutilized by medical physicists and researchers using currently available simulators due to reasons including speed, accuracy, and extensibility constraints. This paper introduces an innovative MRI simulation engine and framework that aims to overcome these issues making available realistic and fast MRIsimulation. Using the CUDA C/C++ programing language, an MRI simulation engine (PhoenixMR), incorporating a Turing-complete virtual machine (VM) to simulate abstract spatiotemporal complexities, was developed. This engine solves a set of time-discrete Bloch equationsusing the symmetric operator splitting technique. An extensible front-end framework package (written in Python) aids the use of PhoenixMR to simplify simulationdevelopment. The PhoenixMR library and front-end codes have been developed and tested. A set of example simulations were performed to demonstrate the ease of use and flexibility of simulation components such as geometrical setup, pulse sequence design, phantom design, and so forth. Initial validation of PhoenixMR is performed by comparing its accuracy and performance against a widely used MRI simulator using identical simulation parameters. Validation results show PhoenixMR simulations are three orders of magnitude faster. There is also strong agreement betweenmodels. A novel MRI simulation platform called PhoenixMR has been introduced. This research tool is designed to be usable by physicists and engineers interested in performing MRI simulations. Examples are shown demonstrating the accuracy, flexibility, and usability of PhoenixMR in several key areas of MRIsimulation.

A Structure-Preserving Semi-implicit IMEX Finite Volume Scheme for Ideal Magnetohydrodynamics at all Mach and Alfvén Numbers

We present a divergence-free semi-implicit finite volume scheme for the simulation of the ideal magnetohydrodynamics (MHD) equations which is stable for large time steps controlled by the local transport speed at all Mach and Alfvén numbers. An operator splitting technique allows to treat the convective terms explicitly while the hydrodynamic pressure and the magnetic field contributions are integrated implicitly, yielding two decoupled linear implicit systems. The linearity of the implicit part is achieved by means of a semi-implicit time linearization. This structure is favorable as second-order accuracy in time can be achieved relying on the class of semi-implicit IMplicit–EXplicit Runge–Kutta (IMEX-RK) methods. In space, implicit cell-centered finite difference operators are designed to discretely preserve the divergence-free property of the magnetic field on three-dimensional Cartesian meshes. The new scheme is also particularly well suited for low Mach number flows and for the incompressible limit of the MHD equations, since no explicit numerical dissipation is added to the implicit contribution and the time step is scale independent. Likewise, highly magnetized flows can benefit from the implicit treatment of the magnetic fluxes, hence improving the computational efficiency of the novel method. The convective terms undergo a shock-capturing second order finite volume discretization to guarantee the effectiveness of the proposed method even for high Mach number flows. The new scheme is benchmarked against a series of test cases for the ideal MHD equations addressing different acoustic and Alfvén Mach number regimes where the performance and the stability of the new scheme is assessed.

Staggered finite volume algorithms for dynamic analysis of saturated soil with low permeability and incompressible pore fluid by using OpenFOAM

Biot dynamic theory forms the basis of dynamic analysis for soil-pore fluid coupling, whose formulations can be simplified by u-p approximation for low frequency phenomena without loss of accuracy. In the past decades, fully-coupled approach and segregated method were two common strategies for Finite Volume (FV) solutions to Biot dynamic theory. However, both methods encounter problems that limit their applicability when soil permeability and pore fluid compressibility reduce to a relatively low level. Under this circumstance, fully-coupled approach requires mixed formulation to satisfy Babuška-Brezzi condition, which complicates numerical implementation. On the other hand, segregated method suffers from divergence and instability, and shows poor performance in computation efficiency. In this article, two staggered FV algorithms incorporating operator splitting technique are devoted to overcoming the aforementioned shortcomings, and they are implemented on the open-source, broadly used and versatile platform for FV calculation, OpenFOAM. Through comparing the computed results of the two algorithms with solutions in the literature, both of the proposed algorithms are shown to reasonably reflect the dynamic responses of porous media. By conducting stability analyses based on the numerical cases, the two algorithms are shown to enlarge the stability region and enable larger time steps for computation.

Simulation of 1D atmospheric pressure dielectric barrier discharge in argon

This work aims at modelling an atmospheric-pressure homogeneous barrier discharge in argon, using a time-dependent 1D fluid model coupled to the electric field and plasmo-chemical kinetic equations. The model is chosen to mimic a discharge when a sinusoidal 1 kV voltage at 10 MHz is applied to the terminals. Energy and mass transfer are considered for a macroscopic fluid representation, while energy transfer in molecular collisions and chemical reactions is treated at the microscopic level. The macroscopic model is represented by a set of coupled partial differential equations. Microscopic effects are studied within a discrete model for electronic and molecular collisions in the frame of ZDPlasKin, a plasma modelling numerical tool. The BOLSIG+ solver is employed in solving the electronic Boltzmann equation. An operator splitting technique is used to separate microscopic and macroscopic models. The spatial and temporal evolution of such species and electron transport parameters are presented and discussed.

Robust dual-angle measurement in magnetization transfer spectroscopy by time-optimal control.

Magnetization transfer spectroscopy relies heavily on the robust determination of relaxation times of nuclei participating in metabolic exchange. Challenges arise due to the use of surface RF coils for transmission (high variation) and the broad resonance band of most X nuclei. These challenges are particularly pronounced when fast mapping methods, such as the dual-angle method, are employed. Consequently, in this work, we develop resonance offset and robust excitation RF pulses for 31P magnetization transfer spectroscopy at 7T through ensemble-based time-optimal control. In our approach, we introduce a cost functional for designing robust pulses, incorporating the full Bloch equations as constraints, which are solved using symmetric operator splitting techniques. The optimal control design of the RF pulses developed demonstrates improved accuracy, desired phase properties, and reduced RF power when applied to dual-angle mapping, thereby improving the precision of exchange-rate measurements, as demonstrated in a preclinical in vivo study quantifying brain creatine kinase activity.

A semi-Lagrangian Splitting framework for the simulation of non-hydrostatic free-surface flows

Incorporation of the free-surface dynamics plays an important role in the numerical modelling of non-hydrostatic flows. In the present work a novel semi-Lagrangian splitting (SLS) scheme for the free-surface Euler system is proposed.Using the operator splitting technique, the three main components emerge: the multilayer Shallow Water Equations (mSWE) (enforcing the free-surface kinematics), the Vertical Remeshing Operator (VRO) (that blends the Eulerian and Lagrangian techniques) and the Pressure Correction Operator (PCO) (that is a classic elliptic equation used in projection-type pressure coupling schemes). The first two are solved by two separate Finite Volume schemes, while the latter is resolved in a Finite Difference manner.In order to assess the SLS’s performance, numerical simulations of periodic waves and solitons propagating over variable bathymetry are performed. The last case concerns the propagation of waves on top of a shear current, showcasing the use of the SLS on rotational flows. The results are compared with experiments and analytic relations and the method is proven to be robust and efficient.

A Multiobjective Optimization Application to Control the Aedes Aegypti Mosquito using a Two-Dimensional Diffusion-Reaction Model

Aedes aegypti is the main vector of multiple diseases, such as dengue, yellow fever, Zika, and chikungunya. Diseases associated with mosquitoes have been growing in recent years, with more than one-third of the world population at risk. Control techniques have been studied to prevent the spread of the Aedes aegypti, such that new formula and how to use adulticides and larvicides, among others. This paper proposes a novel approach in the field of partial differential equations and optimization. We consider a two-dimensional diffusion-reaction model that describes the interaction between aquatic and adult female stages spreading across a domain with parameters dependent on rainfall and temperature. We also formulate a mono-objective and multiobjective optimization problem to minimize the Aedes aegypti populations and the control, considering the application of adulticides and larvicides, using actual data from the city of Lavras/Brazil. The operator splitting technique is used to solve the diffusion-reaction system, coupling finite difference and the fourth-order Runge-Kutta method and optimal solutions were searched by using the Real-Biased Genetic Algorithm and Non-dominated Sorting Genetic Algorithm -II. Numerical results show significant reduction of the Aedes aegypti population.

An interwoven composite tailored finite point method for two dimensional unsteady Burgers' equation

This paper centers on constructing a lucid and utilitarian approach to tackle linear and non-linear two-dimensional partial differential equations. To test the applicability of the proposed algorithm a variant of the classical two-dimensional unsteady Burgers' equation is set up as a testing ground. The method in a nutshell reduces to solving one-dimensional Burgers' equations resulting from the application of appropriate operator splitting techniques in the temporal direction. In solving these one-dimensional Burgers' equations a refined tailored finite point method in conjunction with an apposite linearization to the purpose is employed. The conditional stability, consistency, and convergence of the method are established theoretically and the method is found to be first-order convergent in time and second-order convergent in space. To illustrate the accuracy of the scheme, divers examples have been solved and the results obtained prove that this method is top-notch in terms of cost-cutting and time efficiency through the sufficiency of coarse meshes.

Efficient pricing of options in jump–diffusion models: Novel implicit–explicit methods for numerical valuation

This paper presents novel implicit–explicit Runge–Kutta type methods for numerically simulating partial integro-differential equations that arise when pricing options under jump–diffusion models. These methods offer an alternative approach that avoids the need for numerical or analytical inversion of the coefficient matrix. The pricing of European options is formulated as a partial integro-differential equation, while the pricing of American options are treated as a linear complementarity problem. The developed implicit–explicit Runge–Kutta type method is combined with an operator splitting technique to efficiently solve the linear complementarity problem. Stability and convergence analysis of the proposed methods are established using discrete ℓ2-norm. To validate their efficiency and accuracy, the methods are applied to pricing European and American options under Merton’s and Kou’s models, and the computed results are compared with those reported in the literature.

A numerical approach to the optimal control of thermally convective flows

The optimal control of thermally convective flows is usually modeled by an optimization problem with constraints of Boussinesq equations that consist of the Navier-Stokes equation and an advection-diffusion equation. This optimal control problem is challenging from both theoretical analysis and algorithmic design perspectives. For example, the nonlinearity and coupling of fluid flows and energy transports prevent direct applications of gradient type algorithms in practice. In this paper, we propose an efficient numerical method to solve this problem based on the operator splitting and optimization techniques. In particular, we employ the Marchuk-Yanenko method leveraged by the L2-projection for the time discretization of the Boussinesq equations so that the Boussinesq equations are decomposed into some easier linear equations without any difficulty in deriving the corresponding adjoint system. Consequently, at each iteration, four easy linear advection-diffusion equations and two degenerated Stokes equations at each time step are needed to be solved for computing a gradient. Then, we apply the Bercovier-Pironneau finite element method for space discretization, and design a BFGS type algorithm for solving the fully discretized optimal control problem. We look into the structure of the problem, and design a meticulous strategy to seek step sizes for the BFGS efficiently. Efficiency of the numerical approach is promisingly validated by the results of some preliminary numerical experiments.

Uncoupling Techniques for Multispecies Diffusion–Reaction Model

We consider the multispecies model described by a coupled system of diffusion–reaction equations, where the coupling and nonlinearity are given in the reaction part. We construct a semi-discrete form using a finite volume approximation by space. The fully implicit scheme is used for approximation by time, which leads to solving the coupled nonlinear system of equations at each time step. This paper presents two uncoupling techniques based on the explicit–implicit scheme and the operator-splitting method. In the explicit–implicit scheme, we take the concentration of one species in coupling term from the previous time layer to obtain a linear uncoupled system of equations. The second approach is based on the operator-splitting technique, where we first solve uncoupled equations with the diffusion operator and then solve the equations with the local reaction operator. The stability estimates are derived for both proposed uncoupling schemes. We present a numerical investigation for the uncoupling techniques with varying time step sizes and different scales of the diffusion coefficient.

RAM: Rapid Advection Algorithm on Arbitrary Meshes

The study of many astrophysical flows requires computational algorithms that can capture high Mach number flows, while resolving a large dynamic range in spatial and density scales. In this paper we present a novel method, RAM: Rapid Advection Algorithm on Arbitrary Meshes. RAM is a time-explicit method to solve the advection equation in problems with large bulk velocity on arbitrary computational grids. In comparison with standard upwind algorithms, RAM enables advection with larger time steps and lower truncation errors. Our method is based on the operator splitting technique and conservative interpolation. Depending on the bulk velocity and resolution, RAM can decrease the numerical cost of hydrodynamics by more than one order of magnitude. To quantify the truncation errors and speed-up with RAM, we perform one- and two-dimensional hydrodynamics tests. We find that the order of our method is given by the order of the conservative interpolation and that the effective speed-up is in agreement with the relative increment in time step. RAM will be especially useful for numerical studies of disk−satellite interaction, characterized by high bulk orbital velocities and nontrivial geometries. Our method dramatically lowers the computational cost of simulations that simultaneously resolve the global disk and potential well inside the Hill radius of the secondary companion.

A compact ADI finite difference method for 2D reaction–diffusion equations with variable diffusion coefficients

Reaction–diffusion systems on a spatially heterogeneous domain have been widely used to model various biological applications. However, solving such partial differential equations (PDEs) analytically is rarely possible. Therefore, efficient and accurate numerical methods for solving such PDEs are desired. In this paper, we apply the well-known Padé approximation-based operator splitting techniques and develop a fourth-order compact alternative directional implicit (ADI) scheme. The new scheme is compact and fourth-order accurate in space. Combined with the Richardson extrapolation, the method can be improved to fourth-order accuracy in time. Stability analysis shows that the method is unconditionally stable; thus, a large time step can be used to improve the overall computational efficiency. Numerical examples have also demonstrated the new scheme’s high efficiency and high order accuracy.

An investigation of the multi-mode Richtmyer-Meshkov instability at a gas/HE interface using Pagosa

In this work, we present a hydrocode Pagosa and explore the Richtmyer-Meshkov Instability (RMI) at an air/high explosive (HE) interface for the first time that is important but has not received much attention yet in the high explosive safety field. Thus, the presented Pagosa can be expected to predict the whole deflagration-to-detonation transition (DDT) process in future. In Pagosa, spatial discretization is implemented on cubic staggered grids by computing different variables at the vertex and the cell center, respectively, a special operator-splitting technique is employed to reduce the computational cost, and an artificial viscosity is added to handle the discontinuous shock waves in our simulations. To quantitatively evaluate the capability of Pagosa to solve these kinds of instabilities, the single mode Rayleigh-Taylor instability (RTI) and the multimode RMI at an air/SF 6 interface are first performed, respectively. The Pagosa results are compared with the related numerical solutions in the existing references and the experimental result. Moreover, a theoretical derivation of growth of RTI is also provided based on our numerical method. Subsequently, we explore the multi-mode RMI at an air/HE interface as well as the effects of several factors using Pagosa. Numerical results show that Pagosa is a powerful toolset to generate the right structures and the amplitude of RTI and RMI at an air/SF 6 interface. The solid HE can be penetrated by a strong shock wave and forms RMI deformations. The RMI at an air/HE interface behaves very different than at an air/SF 6 interface, periodic, decreased oscillation is observed due to material character, and is very sensitive to the initial simulation settings, that is, a tiny change in physical quantities will lead to a remarkable RMI structure, which is also observed in a shock bubble interaction [41] . The findings in this work are significant, and will present a new insight for the high explosive field. • An enhanced model Pagosa is presented with some special techniques. • Pagosa's capability to accurately capture the RTI at an air/SF 6 interface. • Pagosa's ability to accurately predict the RMI at an air/SF 6 interface. • The RMI at an air/HE interface is different from that at an air/SF 6 interface. • Effects of several factors on RMI at an air/HE interface.

A DLM/FD method for simulating balls settling in Oldroyd-B viscoelastic fluids

In this article we present a numerical method for simulating the sedimentation of balls in a three-dimensional channel filled with an Oldroyd-B fluid. We have combined a distributed Lagrange multiplier/fictitious domain method with a factorization approach from Lozinski and Owens (2003) [22] via an operator splitting technique. The method is validated by comparing obtained results with those reported in literature. Then the study of fluid elasticity on the formation of ball chain (up to six balls) in Oldroyd-B fluids is presented to demonstrate the capability of our methods. For a higher value of elastic number, a longer chain of balls can be formed while settling in an Oldroyd-B fluid as observed experimentally.

Accelerated Sparse Recovery via Gradient Descent with Nonlinear Conjugate Gradient Momentum

This paper applies an idea of adaptive momentum for the nonlinear conjugate gradient to accelerate optimization problems in sparse recovery. Specifically, we consider two types of minimization problems: a (single) differentiable function and the sum of a non-smooth function and a differentiable function. In the first case, we adopt a fixed step size to avoid the traditional line search and establish the convergence analysis of the proposed algorithm for a quadratic problem. This acceleration is further incorporated with an operator splitting technique to deal with the non-smooth function in the second case. We use the convex $$\ell _1$$ and the nonconvex $$\ell _1-\ell _2$$ functionals as two case studies to demonstrate the efficiency of the proposed approaches over traditional methods.

A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology.

We present a simple, yet efficient adaptive time stepping scheme for cardiac electrophysiology (EP) simulations based on standard operator splitting techniques. The general idea is to exploit the relation between the splitting error and the reaction's magnitude-found in a previous one-dimensional analytical study by Spiteri and Ziaratgahi-to construct the new time step controller for three-dimensional problems. Accordingly, we propose to control the time step length of the operator splitting scheme as a function of the reaction magnitude, in addition to the common approach of adapting the reaction time step. This conforms with observations in numerical experiments supporting the need for a significantly smaller time step length during depolarization than during repolarization. The proposed scheme is compared with classical proportional-integral-differential controllers using state-of-the-art error estimators, which are also presented in details as they have not been previously applied in the context of cardiac EP with operator splitting techniques. Benchmarks show that choosing the time step as a sigmoidal function of the reaction magnitude is highly efficient and full cardiac cycles can be computed with precision even in a realistic biventricular setup. The proposed scheme outperforms common adaptive time stepping techniques, while depending on fewer tuning parameters.