David Pardo

Nathan Collier, Lisandro Dalcin, David Pardo et al.

In this paper we study how the use of a more continuous set of basis functions affects the cost of solving systems of linear equations resulting from a discretized Galerkin weak form. Specifically, we compare performance of linear solvers when discretizing using $C^0$ B-splines, which span traditional finite element spaces, and $C^{p-1}$ B-splines, which represent maximum continuity. We provide theoretical estimates for the increase in cost of the matrix-vector product as well as for the construction and application of black-box preconditioners. We accompany these estimates with numerical results and study their sensitivity to various grid parameters such as element size $h$ and polynomial order of approximation $p$. Finally, we present timing results for a range of preconditioning options for the Laplace problem. We conclude that the matrix-vector product operation is at most $\slfrac{33p^2}{8}$ times more expensive for the more continuous space, although for moderately low $p$, this number is significantly reduced. Moreover, if static condensation is not employed, this number further reduces to at most a value of 8, even for high $p$. Preconditioning options can be up to $p^3$ times more expensive to setup, although this difference significantly decreases for some popular preconditioners such as Incomplete LU factorization.

Carlos Uriarte, David Pardo, Ignacio Muga et al.

Residual minimization is a widely used technique for solving Partial Differential Equations in variational form. It minimizes the dual norm of the residual, which naturally yields a saddle-point (min-max) problem over the so-called trial and test spaces. In the context of neural networks, we can address this min-max approach by employing one network to seek the trial minimum, while another network seeks the test maximizers. However, the resulting method is numerically unstable as we approach the trial solution. To overcome this, we reformulate the residual minimization as an equivalent minimization of a Ritz functional fed by optimal test functions computed from another Ritz functional minimization. We call the resulting scheme the Deep Double Ritz Method (D$^2$RM), which combines two neural networks for approximating trial functions and optimal test functions along a nested double Ritz minimization strategy. Numerical results on different diffusion and convection problems support the robustness of our method, up to the approximation properties of the networks and the training capacity of the optimizers.

Judit Muñoz-Matute, David Pardo, Victor M. Calo et al.

Variational space-time formulations for Partial Differential Equations have been of great interest in the last decades. While it is known that implicit time marching schemes have variational structure, the Galerkin formulation of explicit methods in time remains elusive. In this work, we prove that the explicit Runge-Kutta methods can be expressed as discontinuous Petrov-Galerkin methods both in space and time. We build trial and test spaces for the linear diffusion equation that lead to one, two, and general stage explicit Runge-Kutta methods. This approach enables us to design explicit time-domain (goal-oriented) adaptive algorithms

Maciej Paszynski, Victor Calo, David Pardo

This paper describes a direct solver algorithm for a sequence of finite element meshes that are h-refined towards one or several point singularities. For such a sequence of grids, the solver delivers linear computational cost O(N) in terms of CPU time and memory with respect to the number of unknowns N. The linear computational cost is achieved by utilizing the recursive structure provided by the sequence of h-adaptive grids with a special construction of the elimination tree that allows for reutilization of previously computed partial LU factorizations over the entire unrefined part of the computational mesh. The reutilization technique reduces the computational cost of the entire sequence of h-refined grids from O(N^2) down to O(N). Theoretical estimates are illustrated with numerical results on two- and three-dimensional model problems exhibiting one or several point singularities.

Nathan Collier, David Pardo, Maciej Paszynski et al.

The multi-frontal direct solver is the state-of-the-art algorithm for the direct solution of sparse linear systems. This paper provides computational complexity and memory usage estimates for the application of the multi-frontal direct solver algorithm on linear systems resulting from B-spline-based isogeometric finite elements, where the mesh is a structured grid. Specifically we provide the estimates for systems resulting from $C^{p-1}$ polynomial B-spline spaces and compare them to those obtained using $C^0$ spaces.

Ángel J. Omella, David Pardo

We introduce an $r-$adaptive algorithm to solve Partial Differential Equations using a Deep Neural Network. The proposed method restricts to tensor product meshes and optimizes the boundary node locations in one dimension, from which we build two- or three-dimensional meshes. The method allows the definition of fixed interfaces to design conforming meshes, and enables changes in the topology, i.e., some nodes can jump across fixed interfaces. The method simultaneously optimizes the node locations and the PDE solution values over the resulting mesh. To numerically illustrate the performance of our proposed $r-$adaptive method, we apply it in combination with a collocation method, a Least Squares Method, and a Deep Ritz Method. We focus on the latter to solve one- and two-dimensional problems whose solutions are smooth, singular, and/or exhibit strong gradients.

Tomas Teijeiro, Jamie M. Taylor, Ali Hashemian et al.

We propose the use of machine learning techniques to find optimal quadrature rules for the construction of stiffness and mass matrices in isogeometric analysis (IGA). We initially consider 1D spline spaces of arbitrary degree spanned over uniform and non-uniform knot sequences, and then the generated optimal rules are used for integration over higher-dimensional spaces using tensor product sense. The quadrature rule search is posed as an optimization problem and solved by a machine learning strategy based on gradient-descent. However, since the optimization space is highly non-convex, the success of the search strongly depends on the number of quadrature points and the parameter initialization. Thus, we use a dynamic programming strategy that initializes the parameters from the optimal solution over the spline space with a lower number of knots. With this method, we found optimal quadrature rules for spline spaces when using IGA discretizations with up to 50 uniform elements and polynomial degrees up to 8, showing the generality of the approach in this scenario. For non-uniform partitions, the method also finds an optimal rule in a reasonable number of test cases. We also assess the generated optimal rules in two practical case studies, namely, the eigenvalue problem of the Laplace operator and the eigenfrequency analysis of freeform curved beams, where the latter problem shows the applicability of the method to curved geometries. In particular, the proposed method results in savings with respect to traditional Gaussian integration of up to 44% in 1D, 68% in 2D, and 82% in 3D spaces.

Mohammad Mahdi Abedi, David Pardo, Tariq Alkhalifah

Standard physics-informed neural networks (PINNs) struggle to simulate highly oscillatory Helmholtz solutions in heterogeneous media because pointwise minimization of second-order PDE residuals is computationally expensive, biased toward smooth solutions, and requires artificial absorbing boundary layers to restrict the solution. To overcome these challenges, we introduce a Green-Integral (GI) neural solver for the acoustic Helmholtz equation. It departs from the PDE-residual-based formulation by enforcing wave physics through an integral representation that imposes a nonlocal constraint. Oscillatory behavior and outgoing radiation are encoded directly through the integral kernel, eliminating second-order spatial derivatives and enforcing physical solutions without additional boundary layers. Theoretically, optimizing this GI loss via a neural network acts as a spectrally tuned preconditioned iteration, enabling convergence in heterogeneous media where the classical Born series diverges. By exploiting FFT-based convolution to accelerate the GI loss evaluation, our approach substantially reduces GPU memory usage and training time. However, this efficiency relies on a fixed regular grid, which can limit local resolution. To improve local accuracy in strong scattering regions, we also propose a hybrid GI+PDE loss, enforcing a lightweight Helmholtz residual at a small number of nonuniformly sampled collocation points. We evaluate our method on seismic benchmark models characterized by structural contrasts and subwavelength heterogeneity at frequencies up to 20Hz. GI-based training consistently outperforms PDE-based PINNs, reducing computational cost by over a factor of ten. In models with localized scattering, the hybrid loss yields the most accurate reconstructions, providing a stable, efficient, and physically grounded alternative.

Alejandro Duque, Paulina Sepúlveda, Carlos Uriarte et al.

This work presents a robust, energy-based deep learning framework for solving transmission problems in heterogeneous media, including cases with discontinuous material scenarios. We introduce a weighted First-Order System Least-Squares (FOSLS) formulation involving an energy-norm Poincaré constant and prove its equivalence to a natural energy norm of the underlying equations, with constants independent of material parameters. As a result, the optimization landscape remains aligned with a meaningful error approximation even under high material contrast, where standard neural network losses often deteriorate. We further prove that the FOSLS formulation, together with its integral-loss representation, exhibits a passive variance reduction property, whereby the gradient variance progressively decreases as the loss diminishes, in contrast to methods such as VPINNs and Deep Ritz. From a numerical standpoint, we adopt a reduced-order perspective by constructing a low-dimensional space described by a neural network. The optimal coefficients are computed via a least-squares solver, and the space is subsequently improved through gradient-based updates. By selecting the activation function ReQU, the method mitigates the spurious overshoots typically observed in smooth networks when approximating discontinuities. Numerical experiments in 1D and 2D interface settings corroborate these findings.

Jesus Gonzalez-Sieiro, David Pardo, Vincenzo Nava et al.

We propose a method for reducing the spatial discretization error of coarse computational fluid dynamics (CFD) problems by enhancing the quality of low-resolution simulations using deep learning. We feed the model with fine-grid data after projecting it to the coarse-grid discretization. We substitute the default differencing scheme for the convection term by a feed-forward neural network that interpolates velocities from cell centers to face values to produce velocities that approximate the down-sampled fine-grid data well. The deep learning framework incorporates the open-source CFD code OpenFOAM, resulting in an end-to-end differentiable model. We automatically differentiate the CFD physics using a discrete adjoint code version. We present a fast communication method between TensorFlow (Python) and OpenFOAM (c++) that accelerates the training process. We applied the model to the flow past a square cylinder problem, reducing the error from 120% to 25% in the velocity for simulations inside the training distribution compared to the traditional solver using an x8 coarser mesh. For simulations outside the training distribution, the error reduction in the velocities was about 50%. The training is affordable in terms of time and data samples since the architecture exploits the local features of the physics.

Ibai Ramirez, Joel Pino, David Pardo et al.

Transformers are crucial for reliable and efficient power system operations, particularly in supporting the integration of renewable energy. Effective monitoring of transformer health is critical to maintain grid stability and performance. Thermal insulation ageing is a key transformer failure mode, which is generally tracked by monitoring the hotspot temperature (HST). However, HST measurement is complex, costly, and often estimated from indirect measurements. Existing HST models focus on space-agnostic thermal models, providing worst-case HST estimates. This article introduces a spatio-temporal model for transformer winding temperature and ageing estimation, which leverages physics-based partial differential equations (PDEs) with data-driven Neural Networks (NN) in a Physics Informed Neural Networks (PINNs) configuration to improve prediction accuracy and acquire spatio-temporal resolution. The computational accuracy of the PINN model is improved through the implementation of the Residual-Based Attention (PINN-RBA) scheme that accelerates the PINN model convergence. The PINN-RBA model is benchmarked against self-adaptive attention schemes and classical vanilla PINN configurations. For the first time, PINN based oil temperature predictions are used to estimate spatio-temporal transformer winding temperature values, validated through PDE numerical solution and fiber optic sensor measurements. Furthermore, the spatio-temporal transformer ageing model is inferred, which supports transformer health management decision-making. Results are validated with a distribution transformer operating on a floating photovoltaic power plant.

Mohammad Mahdi Abedi, David Pardo, Tariq Alkhalifah

Physics-Informed Neural Networks (PINNs) have gained increasing attention for solving partial differential equations, including the Helmholtz equation, due to their flexibility and mesh-free formulation. However, their low-frequency bias limits their accuracy and convergence speed for high-frequency wavefield simulations. To alleviate these problems, we propose a simplified PINN framework that incorporates Gabor functions, designed to capture the oscillatory and localized nature of wavefields more effectively. Unlike previous attempts that rely on auxiliary networks to learn Gabor parameters, we redefine the network's task to map input coordinates to a custom Gabor coordinate system, simplifying the training process without increasing the number of trainable parameters compared to a simple PINN. We validate the proposed method across multiple velocity models, including the complex Marmousi and Overthrust models, and demonstrate its superior accuracy, faster convergence, and better robustness features compared to both traditional PINNs and earlier Gabor-based PINNs. Additionally, we propose an efficient integration of a Perfectly Matched Layer (PML) to enhance wavefield behavior near the boundaries. These results suggest that our approach offers an efficient and accurate alternative for scattered wavefield modeling and lays the groundwork for future improvements in PINN-based seismic applications.

Pablo Herrera, Jamie M. Taylor, Carlos Uriarte et al.

Solving Partial Differential Equations (PDEs) using neural networks presents different challenges, including integration errors and spectral bias, often leading to poor approximations. In addition, standard neural network-based methods, such as Physics-Informed Neural Networks (PINNs), often lack stability when dealing with PDEs characterized by low-regularity solutions. To address these limitations, we introduce the Ritz--Uzawa Neural Networks (RUNNs) framework, an iterative methodology to solve strong, weak, and ultra-weak variational formulations. Rewriting the PDE as a sequence of Ritz-type minimization problems within a Uzawa loop provides an iterative framework that, in specific cases, reduces both bias and variance during training. We demonstrate that the strong formulation offers a passive variance reduction mechanism, whereas variance remains persistent in weak and ultra-weak regimes. Furthermore, we address the spectral bias of standard architectures through a data-driven frequency tuning strategy. By initializing a Sinusoidal Fourier Feature Mapping based on the Normalized Cumulative Power Spectral Density (NCPSD) of previous residuals or their proxies, the network dynamically adapts its bandwidth to capture high-frequency components and severe singularities. Numerical experiments demonstrate the robustness of RUNNs, accurately resolving highly oscillatory solutions and successfully recovering a discontinuous $L^2$ solution from a distributional $H^{-2}$ source -- a scenario where standard energy-based methods fail.

Ibai Ramirez, Jokin Alcibar, Joel Pino et al.

Scientific Machine Learning (SciML) integrates physics and data into the learning process, offering improved generalization compared with purely data-driven models. Despite its potential, applications of SciML in prognostics remain limited, partly due to the complexity of incorporating partial differential equations (PDEs) for ageing physics and the scarcity of robust uncertainty quantification methods. This work introduces a Bayesian Physics-Informed Neural Network (B-PINN) framework for probabilistic prognostics estimation. By embedding Bayesian Neural Networks into the PINN architecture, the proposed approach produces principled, uncertainty-aware predictions. The method is applied to a transformer ageing case study, where insulation degradation is primarily driven by thermal stress. The heat diffusion PDE is used as the physical residual, and different prior distributions are investigated to examine their impact on predictive posterior distributions and their ability to encode a priori physical knowledge. The framework is validated against a finite element model developed and tested with real measurements from a solar power plant. Results, benchmarked against a dropout-PINN baseline, show that the proposed B-PINN delivers more reliable prognostic predictions by accurately quantifying predictive uncertainty. This capability is crucial for supporting robust and informed maintenance decision-making in critical power assets.

Mohammad Mahdi Abedi, David Pardo, Tariq Alkhalifah

Physics-Informed Neural Networks (PINNs) have shown promise in solving partial differential equations (PDEs), including the frequency-domain Helmholtz equation. However, standard training of PINNs using gradient descent (GD) suffers from slow convergence and instability, particularly for high-frequency wavefields. For scattered acoustic wavefield simulation based on Helmholtz equation, we derive a hybrid optimization framework that accelerates training convergence by embedding a least-squares (LS) solver directly into the GD loss function. This formulation enables optimal updates for the linear output layer. Our method is applicable with or without perfectly matched layers (PML), and we provide practical tensor-based implementations for both scenarios. Numerical experiments on benchmark velocity models demonstrate that our approach achieves faster convergence, higher accuracy, and improved stability compared to conventional PINN training. In particular, our results show that the LS-enhanced method converges rapidly even in cases where standard GD-based training fails. The LS solver operates on a small normal matrix, ensuring minimal computational overhead and making the method scalable for large-scale wavefield simulations.

Kyubo Noh, David Pardo, Carlos Torres-Verdin

Deep Learning (DL) inversion is a promising method for real time interpretation of logging while drilling (LWD) resistivity measurements for well navigation applications. In this context, measurement noise may significantly affect inversion results. Existing publications examining the effects of measurement noise on DL inversion results are scarce. We develop a method to generate training data sets and construct DL architectures that enhance the robustness of DL inversion methods in the presence of noisy LWD resistivity measurements. We use two synthetic resistivity models to test three approaches that explicitly consider the presence of noise: (1) adding noise to the measurements in the training set, (2) augmenting the training set by replicating it and adding varying noise realizations, and (3) adding a noise layer in the DL architecture. Numerical results confirm that the three approaches produce a denoising effect, yielding better inversion results in both predicted earth model and measurements compared not only to the basic DL inversion but also to traditional gradient based inversion results. A combination of the second and third approaches delivers the best results. The proposed methods can be readily generalized to multi dimensional DL inversion.

Sergey Alyaev, Mostafa Shahriari, David Pardo et al.

Modern geosteering is heavily dependent on real-time interpretation of deep electromagnetic (EM) measurements. We present a methodology to construct a deep neural network (DNN) model trained to reproduce a full set of extra-deep EM logs consisting of 22 measurements per logging position. The model is trained in a 1D layered environment consisting of up to seven layers with different resistivity values. A commercial simulator provided by a tool vendor is used to generate a training dataset. The dataset size is limited because the simulator provided by the vendor is optimized for sequential execution. Therefore, we design a training dataset that embraces the geological rules and geosteering specifics supported by the forward model. We use this dataset to produce an EM simulator based on a DNN without access to the proprietary information about the EM tool configuration or the original simulator source code. Despite employing a relatively small training set size, the resulting DNN forward model is quite accurate for the considered examples: a multi-layer synthetic case and a section of a published historical operation from the Goliat Field. The observed average evaluation time of 0.15 ms per logging position makes it also suitable for future use as part of evaluation-hungry statistical and/or Monte-Carlo inversion algorithms within geosteering workflows.