Deep learning-based hybrid iterative methods (DL-HIM) have emerged as a promising approach for designing fast neural solvers to tackle large-scale sparse linear systems. DL-HIM combine the smoothing effect of simple iterative methods with the spectral bias of neural networks, which allows them to effectively eliminate both high-frequency and low-frequency error components. However, their efficiency may decrease if simple iterative methods can not provide effective smoothing, making it difficult for the neural network to learn mid-frequency and high-frequency components. This paper first conducts a convergence analysis for general DL-HIM from a spectral viewpoint, concluding that under reasonable assumptions, DL-HIM exhibit a convergence rate independent of grid size $h$ and physical parameters $\boldsymbol{\mu}$. To meet these assumptions, we design a neural network from an eigen perspective, focusing on learning the eigenvalues and eigenvectors corresponding to error components that simple iterative methods struggle to eliminate. Specifically, the eigenvalues are learned by a meta subnet, while the eigenvectors are approximated using Fourier modes with a transition matrix provided by another meta subnet. The resulting DL-HIM, termed the Fourier Neural Solver (FNS), can be trained to achieve a convergence rate independent of PDE parameters and grid size within a local neighborhood of the training scale by designing a loss function that ensures the neural network complements the smoothing effect of the damped Jacobi iterative methods. We verify the performance of FNS on five types of linear parametric PDEs.

This work considers to numerically solve a subdiffusion equation involving constant time delay $\tau$ and Riemann-Liouville fractional derivative. First, a fully discrete finite element scheme is developed for the considered problem under the symmetric graded time mesh, where the Caputo fractional derivative is approximated via the L1 formula, while the Riemann-Liouville integral is discretized using the fractional right rectangular rule. Under the assumption that the exact solution has low regularities at $t=0$ and $\tau$, the local truncation errors of both the L1 formula and the fractional right rectangular rule are analyzed. It is worth noting that, by setting the mesh parameter $r=1$, the symmetric graded time mesh will degenerate to a uniform mesh. Consequently, we proceed to discuss the stability and convergence of the proposed numerical scheme under two scenarios. For the uniform time mesh, by introducing a discrete sequence $\{P_k\}$, the unconditional stability and local time error estimate for the developed scheme is established. Conversely, on the symmetric graded time mesh, through the introduction of a discrete fractional Gronwall inequality, the stability and globally optimal time error estimate can be obtained. Finally, some numerical tests are presented to validate the theoretical results.

The phase field model is a widely used mathematical approach for describing crack propagation in continuum damage fractures. In the context of phase field fracture simulations, adaptive finite element methods (AFEM) are often employed to address the mesh size dependency of the model. However, existing AFEM approaches for this application frequently rely on heuristic adjustments and empirical parameters for mesh refinement. In this paper, we introduce an adaptive finite element method based on a recovery type posteriori error estimates approach grounded in theoretical analysis. This method transforms the gradient of the numerical solution into a smoother function space, using the difference between the recovered gradient and the original numerical gradient as an error indicator for adaptive mesh refinement. This enables the automatic capture of crack propagation directions without the need for empirical parameters. We have implemented this adaptive method for the Hybrid formulation of the phase field model using the open-source software package FEALPy. The accuracy and efficiency of the proposed approach are demonstrated through simulations of classical 2D and 3D brittle fracture examples, validating the robustness and effectiveness of our implementation.

This paper proposes an efficient general alternating-direction implicit (GADI) framework for solving large sparse linear systems. The convergence property of the GADI framework is discussed. Most of the existing ADI methods can be viewed as particular schemes of the developed framework. Meanwhile the GADI framework can derive new ADI methods. Moreover, as the algorithm efficiency is sensitive to the splitting parameters, we offer a data-driven approach, the Gaussian process regression (GPR) method based on the Bayesian inference, to predict the GADI framework's relatively optimal parameters. The GPR method requires a small training data set to learn the regression prediction mapping, which has sufficient accuracy and high generalization capability. It allows us to efficiently solve linear systems with a one-shot computation, and does not require any repeated computations to obtain relatively optimal splitting parameters. Finally, we use the three-dimensional convection-diffusion equation and continuous Sylvester matrix equation to examine the performance of our proposed methods. Numerical results demonstrate that the proposed framework is faster tens to thousands of times than the existing ADI methods, such as (inexact) Hermitian and skew-Hermitian splitting type methods in which the consumption of obtaining relatively optimal splitting parameters is ignored. Due to the efficiency of the developed methods, we can solve much larger linear systems which these existing ADI methods have been not reached.

Block copolymers provide a wonderful platform in studying the soft condensed matter systems. Many fascinating ordered structures have been discovered in bulk and confined systems. Among various theories, the self-consistent field theory (SCFT) has been proven to be a powerful tool for studying the equilibrium ordered structures. Many numerical methods have been developed to solve the SCFT model. However, most of these focus on the bulk systems, and little work on the confined systems, especially on general curved surfaces. In this work, we developed a linear surface finite element method, which has a rigorous mathematical theory to guarantee numerical precsion, to study the self-assembled phases of block copolymers on general curved surfaces based on the SCFT. Furthermore, to capture the consistent surface for a given self-assembled pattern, an adaptive approach to optimize the size of the general curved surface has been proposed. To demonstrate the power of this approach, we investigate the self-assembled patterns of diblock copolymers on several distinct curved surfaces, including five closed surfaces and an unclosed surface. Numerical results illustrate the efficiency of the proposed method. The obtained ordered structures are consistent with the previous results on standard surfaces, such as sphere and torus. Certainly, the proposed numerical framework has the capability of studying the phase behaviors on general surfaces precisely.

This paper aims to develop an efficient adaptive finite element method for the second-order elliptic problem. Although the theory for adaptive finite element methods based on residual-type a posteriori error estimator and bisection refinement has been well established, in practical computations, the use of non-asymptotic exact of error estimator and the excessive number of adaptive iteration steps often lead to inefficiency of the adaptive algorithm. We propose an efficient adaptive finite element method based on high-accuracy techniques including the superconvergence recovery technique and high-quality mesh optimization. The centroidal Voronoi Delaunay triangulation mesh optimization is embedded in the mesh adaption to provide high-quality mesh, and then assure that the superconvergence property of the recovered gradient and the asymptotical exactness of the error estimator. A tailored adaptive strategy, which could generate high-quality meshes with a target number of vertices, is developed to ensure the adaptive computation process terminated within $7$ steps. The effectiveness and robustness of the adaptive algorithm is numerically demonstrated.

Building upon the recent work of Teso and Plociniczak (2025) regarding L1 discretization errors for the Caputo derivative in H\"{o}lder spaces, this study extends the analysis to higher-order discretization errors within the same functional framework. We first investigate truncation errors for the L2 and L1-2 methods, which approximate the Caputo derivative via piecewise quadratic interpolation. Then we generalize the results to arbitrary high-order discretization. Theoretical analyses reveal a unified error structure across all schemes: the convergence order equals the difference between the smoothness degree of the function space and the fractional derivative order, i.e., order of error = degree of smoothness - order of the derivative. Numerical experiments validate these theoretical findings.

In recent years, topology optimization (TO) has gained widespread attention as a powerful structural design method. However, its application remains challenging due to the deep expertise and extensive development effort required. Traditional TO methods, tightly coupled with computational mechanics like finite element method (FEM), result in intrusive algorithms demanding a comprehensive system understanding. This paper presents SOPTX, a TO package based on FEALPy, which implements a modular architecture that decouples analysis from optimization, supports multiple computational backends (NumPy, PyTorch, JAX), and achieves a non-intrusive design paradigm. Core innovations include: (1) cross-platform design that supports multiple computational backends, enabling efficient algorithm execution on central processing units (CPUs) and flexible acceleration using graphics processing units (GPUs), while leveraging automatic differentiation (AD) technology for efficient sensitivity computation of objective and constraint functions; (2) fast matrix assembly techniques that overcome the performance bottlenecks of traditional numerical integration methods, significantly accelerating finite element computations and enhancing overall efficiency; (3) a modular framework supporting TO problems for arbitrary dimensions and meshes, allowing flexible configuration and extensibility of optimization workflows through a rich library of composable components. Using the density-based method for the classic compliance minimization problem with volume constraints as an example, numerical experiments demonstrate SOPTX's high efficiency in computational speed and memory usage, while showcasing its strong potential for research and engineering applications.

The high Weissenberg number problem has been a persistent challenge in the numerical simulation of viscoelastic fluid flows. This paper presents an improved lattice Boltzmann method for solving viscoelastic flow problems at high Weissenberg numbers. The proposed approach employs two independent two-relaxation-time regularized lattice Boltzmann models to solve the hydrodynamic field and conformation tensor field of viscoelastic fluid flows, respectively. The viscoelastic stress computed from the conformation tensor is directly embedded into the hydrodynamic field using a newly proposed local velocity discretization scheme, thereby avoiding spatial gradient calculations. The constitutive equations are treated as convection-diffusion equations and solved using an improved convection-diffusion model specifically designed for this purpose, incorporating a novel auxiliary source term that eliminates the need for spatial and temporal derivative computations. Additionally, a conservative non-equilibrium bounce-back (CNEBB) scheme is proposed for implementing solid wall boundary conditions in the constitutive equations. The robustness of the present algorithm is validated through a series of benchmark problems. The simplified four-roll mill problem demonstrates that the method effectively improves numerical accuracy and stability in bulk regions containing stress singularities. The Poiseuille flow problem validates the accuracy of the current algorithm with the CNEBB boundary scheme at extremely high Weissenberg numbers (tested up to Wi = 10,000). The flow past a circular cylinder problem confirms the superior stability and applicability of the algorithm for complex curved boundary problems compared to other existing common schemes.

The solution of nonsymmetric but positive definite (NSPD) systems arising from advection-diffusion problems is an important research topic in science and engineering. Balancing domain decomposition by constraints with an adaptive coarse space(adaptive BDDC) constitute a significant class of nonoverlapping domain decomposition methods, commonly used for symmetric positive definite problems. In this paper, we propose an adaptive BDDC method that incorporates a class of edge generalized eigenvalue problems based on prior selected primal constraints to solve NSPD systems from advection-diffusion problems. Compared with the conventional adaptive BDDC method for such systems, the proposed approach further reduces the number of primal unknowns. Numerical experiments show that although the iteration count increases slightly, the overall computational time is significantly reduced.

In this paper, we construct and analyze an energy stable scheme by combining the latest developed scalar auxiliary variable (SAV) approach and linear finite element method (FEM) for phase field crystal (PFC) model, and show rigorously that the scheme is first-order in time and second-order in space for the L 2 and H -1 gradient flow equations. To reduce efficiently computational cost and capture accurately the phase interface, we give a simple adaptive strategy, equipped with a posteriori gradient estimator, i.e. L 2 norm of the recovered gradient. Extensive numerical experiments are presented to verify our theoretical results and to demonstrate the effectiveness and accuracy of our proposed method.

In this work, we design a multi-category inverse design neural network to map ordered periodic structure to physical parameters. The neural network model consists of two parts, a classifier and Structure-Parameter-Mapping (SPM) subnets. The classifier is used to identify structure, and the SPM subnets are used to predict physical parameters for desired structures. We also present an extensible reciprocal-space data augmentation method to guarantee the rotation and translation invariant of periodic structures. We apply the proposed network model and data augmentation method to two-dimensional diblock copolymers based on the Landau-Brazovskii model. Results show that the multi-category inverse design neural network is high accuracy in predicting physical parameters for desired structures. Moreover, the idea of multi-categorization can also be extended to other inverse design problems.

Proximal gradient-based optimization is one of the most common strategies to solve inverse problem of images, and it is easy to implement. However, these techniques often generate heavy artifacts in image reconstruction. One of the most popular refinement methods is to fine-tune the regularization parameter to alleviate such artifacts, but it may not always be sufficient or applicable due to increased computational costs. In this work, we propose a deep geometric incremental learning framework based on the second Nesterov proximal gradient optimization. The proposed end-to-end network not only has the powerful learning ability for high-/low-frequency image features, but also can theoretically guarantee that geometric texture details will be reconstructed from preliminary linear reconstruction. Furthermore, it can avoid the risk of intermediate reconstruction results falling outside the geometric decomposition domains and achieve fast convergence. Our reconstruction framework is decomposed into four modules including general linear reconstruction, cascade geometric incremental restoration, Nesterov acceleration, and post-processing. In the image restoration step, a cascade geometric incremental learning module is designed to compensate for missing texture information from different geometric spectral decomposition domains. Inspired by the overlap-tile strategy, we also develop a post-processing module to remove the block effect in patch-wise-based natural image reconstruction. All parameters in the proposed model are learnable, an adaptive initialization technique of physical parameters is also employed to make model flexibility and ensure converging smoothly. We compare the reconstruction performance of the proposed method with existing state-of-the-art methods to demonstrate its superiority. Our source codes are available at this https URL

Finding index-1 saddle points is crucial for understanding phase transitions. In this work, we propose a simple yet efficient approach, the spring pair method (SPM), to accurately locate saddle points. Without requiring Hessian information, SPM evolves a single pair of spring-coupled particles on the potential energy surface. By cleverly designing complementary drifting and climbing dynamics based on gradient decomposition, the spring pair converges onto the minimum energy path (MEP) and spontaneously aligns its orientation with the MEP tangent, providing a reliable ascent direction for efficient convergence to saddle points. SPM fundamentally differs from traditional surface walking methods, which rely on the eigenvectors of Hessian that may deviate from the MEP tangent, potentially leading to convergence failure or undesired saddle points. The efficiency of SPM for finding saddle points is verified by ample examples, including high-dimensional Lennard-Jones cluster rearrangement and the Landau energy functional involving quasicrystal phase transitions.