The Hadamard product of tensor train (TT) tensors is a fundamental nonlinear operation in scientific computing and data analysis. However, due to its tendency to significantly increase TT ranks, the Hadamard product poses a major computational challenge in TT-tensor-based algorithms. To address this, it is crucial to develop recompression algorithms that mitigate the effects of this rank increase. Existing recompression algorithms require an explicit representation of the Hadamard product, resulting in high computational and storage costs. In this work, we propose a Hadamard-avoiding TT recompression (HaTT) algorithm, which reduces both computational complexity and storage requirements. By leveraging the structure of the Hadamard product of TT tensors and exploiting its Hadamard product-free property, the HaTT algorithm achieves significantly lower complexity compared to existing TT recompression methods. This is confirmed through both complexity analysis and numerical experiments. Furthermore, the HaTT algorithm is applied to solve the Allen-Cahn equation, achieving substantial speedup over existing TT recompression algorithms without sacrificing accuracy.
Abstract. We present a conservation error analysis of non-conservative modified ghost fluid methods (MGFM) for axisymmetric compressible multi-medium flows within a discontinuous Galerkin framework. We propose a corresponding correction strategy. While multi-medium Riemann-solver-based MGFM-type methods (e.g., MGFM, MGFM) excel at simulating strong shock-interface interactions, their inherent non-conservatism can induce mass loss and pressure misalignment in long-time, high-energy-density computations. To address this issue, we theoretically analyze the conservation error introduced by the interface treatment, the cell-wise conservation error and the flux mismatch error. We prove that for MGFM-type methods, the global conservation error CE(M,N) is O(Δt). These errors are then conservatively redistributed to the grid cells adjacent to the interface, thereby eliminating the numerical losses of mass, momentum, and energy. The cell-wise error is redistributed using a local smoothing technique, while the flux mismatch error is allocated according to the wave propagation direction determined by the multi-material Riemann problem. The resulting algorithm integrates high-order discontinuous Galerkin discretizations with one-sided Riemann problems for the symmetry axis, forming a robust numerical framework for axisymmetric multi-medium flows in high-energy-density environments. Numerical experiments demonstrate that the proposed method not only achieves global conservation to machine precision but also effectively eliminates the pressure misalignment at the interface. In long-term, large-scale multi-medium simulations, the predicted key physical quantities are in excellent agreement with experimental data and reference solutions.
In high-temperature plasma physics, strong magnetic fields are essential for confining charged particles. Consequently, classical mathematical models of such systems must take into account the external magnetic field effects. A key governing equation is the magnetized Vlasov–Poisson system, which exhibits multiscale dynamics and rich physical properties. Thus, developing the structure-preserving numerical methods and studying its capability of maintaining these intrinsic properties over long-time simulations is therefore critically important. This paper presents a general framework for constructing and analyzing structure-preserving methods in orthogonal curvilinear coordinates. We prove that in these coordinates the Poisson-bracket structure is retained under appropriate finite element discretizations. However, the resulting Hamiltonian systems in transformed coordinates typically can’t be decomposed as several subsystems which can be solved exactly. To address this, we propose a semi-implicit numerical scheme that can still maintain the favorable stability properties inherited by the system. The effectiveness of the new derived numerical methods is demonstrated in application to strongly magnetized systems, and the rigorous asymptotic stability analysis are provided.
This paper presents a fourth-order Cartesian grid based kernel-free boundary integral (KFBI) method for two and three dimensional boundary value problems and interface problems of general variable coefficients elliptic equation. The method reformulates the two elliptic problems into the corresponding boundary integral equations (BIEs) according to the potential theory, and reinterprets the integrals involved as solutions to some simple interface problems in an extended regular domain. The BIEs are iteratively solved by the GMRES method. During each iteration, the equivalent simple interface problems are solved by a fourth-order Cartesian grid based finite difference method, where the sharp interface conditions are accurately captured through a Taylor expansion based correction technique. The essence of "kernel-free" lies in that the Green functions used here are especially defined and never formulated or computed in the whole calculation process. Numerical experiments for two kinds of elliptic problems in both two and three dimensions are shown to demonstrate the algorithm efficiency and numerical accuracy.
We present a comprehensive review of finite element methods (FEM) for solving fluid–structure interaction (FSI) problems. Our aim is to systematically categorise and compare these methods based on three key aspects: mesh type, coupling strategy, and the choice of solved variables. For each method, we highlight its respective advantages and limitations. In addition, we formulate and compare the finite element weak forms of the various approaches, employing Newton’s method to linearise nonlinear terms at the continuous functional level. We implement six FSI methods in the open-source software package FreeFEM++, and compare these methods based on three benchmarks.
In this paper, we propose two novel explicit uniformly accurate exponential integrators for solving the relativistic charged-particle dynamics under a strong magnetic field, which is commonly utilized in the guiding center theory. The solutions exhibit highly oscillatory behavior in time, characterized by a small parameter 0<ε≪1. For constant magnetic field intensity, we introduce a transformation to filter out the oscillatory term entirely. Combining the two-scale formulation, two explicit exponential integrators with second-order and fourth-order uniform accuracy are constructed and a rigorous convergence analysis is provided within the 4D motion equation framework. The error estimate demonstrates that the accuracy and computational cost of the algorithms are both independent of the magnetic field strength. By re-scaling the time, the methods are extended to handle the general strong magnetic field with varying intensity and direction. Moreover, the proposed numerical schemes are also applicable to scenarios involving the maximal ordering scaling magnetic field. Numerical experiments are carried out to validate the efficiency of the methods across diverse magnetic fields.
Thin liquid films with contact lines are common in nature and in engineering applications. The existence of free boundaries and singularities poses significant challenges for both modeling and computation. In this work, we present a unified framework for the derived and numerical approximation of a fourth-order thin film equation with macroscopic dynamic boundary conditions. The reduced model is systematically derived using the Onsager variational principle in conjunction with lubrication theory, yielding a thermodynamically consistent formulation that accounts for capillarity, gravity, and external force. To solve the resulting free boundary problem, we develop adaptive moving mesh methods based on a discrete Onsager variational principle, including a stabilized semi-implicit scheme to improve computational efficiency. Numerical results confirm the optimal convergence of the proposed methods and accurately capture key wetting behaviors, including contact angle hysteresis on rough substrates. This work offers a robust framework for simulating thin film flows with complex geometries.
The ground states of rotating Bose-Einstein condensates (BECs) can be accurately computed as steady-state solutions of the normalized gradient flow with Lagrange multiplier (GFLM) and the reconfigured continuous normalized gradient flow (RCNGF). Although high-order numerical schemes can reduce the number of iterations to reach steady state, their practical efficiency is often limited by a significantly increased computational cost per iteration. In this paper, we propose a class of exponential time differencing Runge-Kutta schemes that preserve the steady-state solutions of both the GFLM and RCNGF. These schemes employ explicit updates while maintaining low per-iteration computational costs. Our newly developed algorithms outperform conventional backward-forward Euler (BF) schemes in computing the ground states of rotating BECs. Numerical experiments demonstrate that, under the same time steps and convergence criteria, the second-order schemes require only about half the iterations of the BF scheme while achieving higher accuracy. Specifically, the proposed schemes yield lower ground-state energies and smaller residuals in the Euler-Lagrange equations for the eigenvalue problem.
Stabilization techniques offer significant advantages in developing highly stable algorithms for gradient flows; however, the considerable lagging effect arising from different integrations of stabilization terms continues to pose a major challenge. For a class of L² gradient flows with derivative-independent nonlinearity, we propose a unified framework to analyze the rescaled time step and the associated enlargement of the time-step constraint in a class of stabilized single-step schemes. As examples, we consider the implicit-explicit Runge–Kutta schemes, and exponential-time-differencing Runge–Kutta schemes. A unified matrix-vector framework is developed to establish their energy stability. We first demonstrate that the unstabilized scheme can maintain the energy dissipation under a specific time-step constraint. In contrast, the stabilized scheme exhibits energy dissipation for any time step, provided that the stabilization parameter is sufficiently large. By reformulating these stabilized single-step schemes into a class of explicit Runge–Kutta integrators, we characterize the time delay introduced by stabilization and eliminate the lagging phenomenon through a relaxation technique. The maximum rescaled time step shows a considerable enlargement of the time-step constraint compared to the unstabilized scheme. Numerical experiments confirm the delay-free behavior, efficiency, and energy dissipation property of proposed schemes.
Liquid crystal hydrodynamics exhibits complex multiscale phenomena governed by strong nonlinearities and anisotropy, posing significant challenges for computational simulation. In this paper, we develop a series of high-order linear and decoupled implicit-explicit (IMEX) schemes for the Ericksen–Leslie model of nematic liquid crystal flows based on the sphere projection. Only two elliptic equations with constant coefficients at each time step need to be solved, which makes the schemes very efficient in simulating such complex fluid systems. A distinctive feature of the proposed IMEX schemes is their ability to assure unconditional energy stability and the preservation of the length of the director field, pivotal for maintaining the physical fidelity of the simulations. We prove that the numerical solutions remain uniformly bounded without constraints on the time step size. Furthermore, we establish rigorous error estimates with orders ranging from 1 to 5 within a unified framework. Several numerical benchmark simulations, including the three-dimensional case that has not been reported in previous studies, are presented to validate the theoretical results of the proposed schemes and to illustrate the defect dynamics of liquid crystal motion.
The mathematical model of electroporoelasticity is comprised of Maxwell's equations coupled to Biot's equations. Numerical approximations may exhibit Poisson locking when the Lamé constant is large and the exact solution's divergence is small. In this paper, we propose a locking-free finite element method for approximating solutions, derive error estimates and show first-order convergence (which are independent of the Lamé constant). It is shown that the proposed numerical method inherits its asymptotic behavior from the original equations, ensuring that the algorithm is locking-free. We present numerical experiments that illustrate our theoretical findings.
In this paper, based on Rosenbrock methods and relaxation techniques, we introduce a class of implicit-explicit methods called the relaxed implicit-explicit Runge-Kutta-Rosenbrock (RIMEXRKR) methods, which are designed to solve problems with stiff and nonstiff terms. The main advantage of our implicit-explicit methods is that they are monotonicity-preserving/conservative and easy to implement. To show the superior properties, we construct novel two-stage second-order L-stable numerical schemes, and analyze their stability regions. In contrast, most existing implicit-explicit methods achieve the same second-order convergence but require three stages. Finally, through several numerical examples of typical partial differential equations, the numerical results validate the theoretical results.
Computing equilibrium shapes of crystals (ESC) is an intriguing problem arising from materials science. In this paper, we consider an unsupervised deep learning approach by proposing a novel symmetrized deep neural network (SDNN) to address the problem. Within the phase-field framework, this problem can be formulated as minimizing an orientation-dependent interfacial free energy functional, subject to the constraint of a given constant area in 2D or volume in 3D. Subsequently, we utilize SDNN to solve the constrained minimization problem. Numerical findings illustrate that the SDNN can precisely determine the equilibrium interface through the zero level-set of the phase function. Moreover, we incorporate more sophisticated strategies, such as the pre-training, L-BFGS optimizer, and adaptive sampling, to enhance both the accuracy and efficiency of the SDNN. Ample numerical results demonstrate the superiority of the SDNN approach in solving the ESC problem. As an initial step, we are convinced that the proposed SDNN approach holds significant potential for exploring the mystery of the ESC.
In this paper, we propose an efficient numerical algorithm for solving Euler's elastica-based inpainting and segmentation models. These models involve a complex curvature term that is nonsmooth and nonconvex, presenting significant numerical challenges. To address these challenges, we extend previous works to develop an accelerated operator splitting (AOS) algorithm by integrating the inertial extrapolation technique into the operator-splitting method based on Lie scheme and Marchuk-Yanenko discretization. By employing a reweighting technique, the proposed AOS algorithm is easy to be implemented, as it mainly requires soft shrinkage and FFT calculation at each iteration. Numerical experiments on image inpainting and segmentation demonstrate the effectiveness and efficiency of the proposed method, especially the remarkable superiority in both iteration number and running time.
This paper aims to develop the novel well-balanced fifth-order moment-based Hermite weighted essentially non-oscillatory (HWENO) schemes for one- and two-dimensional shallow water equations with non-flat bottom topography. A new CST (constant subtraction technique) pre-balanced form [Yang et al., J. Sci. Comput., 63(2015):678-698] is adopted to design the well-balanced scheme easily. Both the flux gradient and source term in the new CST pre-balanced form approach zero for the steady-state solution. The fifth-order moment-based HWENO approximations, which combine a fourth-degree polynomial with a set of second-degree polynomials convexly in $\mathbb{P}^k$ polynomial space are employed in the spatial reconstruction. The corresponding linear weights can be manually chosen as positive numbers, as long as their sum equals one. It offers a better resolution for the perturbation of steady-state flows. Extensive one- and two-dimensional numerical examples are presented to illustrate the fifth-order accuracy, non-oscillatory, good resolution, and well-balanced properties of the proposed method.
Deep learning has been widely used for solving stochastic partial differential equations, especially for high-dimensional problems. Most existing deep neural networks often adopt data-driven or physics-driven methods for training. Pure data-driven methods easily cause model overfitting, and the physics-driven method may fail to optimize for sophisticated PDE systems. The networks are mostly fully connected layers or deep residual networks, which results in a large number of parameters and poor performance for intricate PDEs. In this paper, we design a novel deep neural network surrogate model for high-dimensional stochastic PDEs which is called DenseLeNet. We combine both data-driven and physical constraints to train our model. This hybrid approach enables us to overcome the limitations of purely data-driven or physics-based methods. We test our method through a stochastic elliptic partial differential equation (SPDE) with different high-dimensional uncertain inputs. Numerical results demonstrate that our approach is more generalizable than current DNN models such as ResNet, DenseNet, as well as GoogLeNet.
Deep learning has gained significant development in the field of scientific computing, especially in its application to solve problems related to differential operators using deep neural networks. However, the utilization of neural networks to solve problems involving singularities still faces challenges. In this paper, we discussed the failure of deep learning methods for the singular variational problems exhibiting the Lavrentiev phenomenon. For such problems, we show that the standard deep Ritz method and some variants fail to detect the singular minimizers. We then introduce a guiding term that renders the neural network to explore solutions as desired during training. Numerical experiments demonstrate that the method achieves much better approximations than the previous methods. Furthermore, we apply the same algorithm to solve problems with regular solutions to show the robustness of the proposed method.
This paper introduces a fully discrete finite element numerical scheme specifically tailored for a ternary time-dependent Ginzburg-Landau-de Gennes mesoscopic model, and provides a thorough analysis. The scheme is constructed using the mass-lumped finite element spatial discretization and the second-order Backward Differentiation Formula (BDF2) temporal discretization, ensuring a robust framework for simulation. To effectively handle the highly complicated energy functional, we employ the convex-concave decomposition technique. Meanwhile, an addition of a Douglas-Dupont regularization term allows one to establish a modified energy stability estimate, which is crucial for the theoretical foundation. The unique solvability and positivity-preserving properties of the proposed scheme are rigorously justified, reinforcing its theoretical validity. To validate these theoretical insights and demonstrate the practical applications, a series of numerical experiments are performed, highlighting the evolution of the MMC hydrogel. Notably, the study also delves into the impact of the stability factor on both energy stability and numerical error, providing a comprehensive evaluation of the numerical performance and reliability.
In this work, we study the behaviors of electromagnetic waves in the one-dimensional Cole-Cole dispersive medium surrounded by perfectly matched layers (C-CPML). A combined scheme with finite difference method in time and discontinuous Galerkin (DG) method in space is developed to discretize the C-CPML model. We carefully design the energy functionals in continuous, semi-discrete, and fully discrete forms, and rigorously establish the energy dissipation property of the corresponding systems. Optimal convergence rates of the semi-discrete and fully discrete systems are established theoretically and verified numerically. In addition to convergence order tests, numerical experiments are conducted to evaluate the performance of the PML in absorbing unphysical wave reflections. Numerical results demonstrate that the PML effectively dissipates waves entering the layer and leads to negligible reflections.
Whole-body hemodynamic simulations provide quantitative tools for patient-specific treatment planning but require scalable computational methods to model complex physiological systems. We develop a parallel solver within the Newton-Krylov-Schwarz framework to simulate 3D blood flow in a patient-specific arterial network with over 500 branches, governed by the incompressible Navier-Stokes equations on unstructured tetrahedral meshes. Our novel Hybrid Parallel Block Iterative Solver (HPBIS) and its Reduced and Reused variant (RR-HPBIS), optimized for the Sunway TaihuLight’s many-core architecture, employ batched data transfers and computation-communication overlap. Numerical experiments show RR-HPBIS achieves a 25.2x speedup over a homogeneous baseline and scales to 16,384 CPU cores for a 197-million-tetrahedron mesh, yielding a 5.1x speedup at 32% efficiency. The solver captures systemic hemodynamics, including spiral flow and wall shear stress, critical for vascular physiology. Despite simplified boundary conditions, this method advances high-resolution computational modeling for clinical applications.
