ConDiff: A Challenging Dataset for Neural Solvers of Partial Differential Equations

## Abstract

Abstract We present ConDiff, a novel dataset for scientific machine learning. ConDiff focuses on the diffusion equation with varying coefficients, a fundamental problem in many applications of parametric partial differential equations (PDEs). The main novelty of the proposed dataset is that we consider discontinuous coefficients with high contrast. These coefficient functions are sampled from a selected set of distributions. This class of problems is not only of great academic interest, but is also the basis for describing various environmental and industrial problems. In this way, ConDiff shortens the gap with real-world problems while remaining fully synthetic and easy to use. ConDiff consists of a diverse set of diffusion equations with coefficients covering a wide range of contrast levels and heterogeneity with a measurable complexity metric for clearer comparison between different coefficient functions. We baseline ConDiff on standard deep learning models in the field of scientific machine learning. By providing a large number of problem instances, each with its own coefficient function and right-hand side, we hope to encourage the development of novel physics-based deep learning approaches, such as neural operators and physics-informed neural networks, ultimately driving progress towards more accurate and efficient solutions of complex PDE problems.

## 1 Introduction

In recent years, machine learning techniques have emerged as a promising approach to solving PDEs, offering a new perspective in scientific computing. Machine learning algorithms, especially those based on neural networks, have demonstrated success in approximating complex functions and physical phenomena. Neural networks can provide more efficient and scalable methods compared to traditional numerical methods, which can be computationally expensive and limited by the dimensionality of the problem to be solved. Approaches using physical losses Karniadakis et al. (), operator learning Li et al. (), symmetries incorporation Wang et al. (), data-driven discretization Bar-Sinai et al. () lead to more physically meaningful solutions and gave neural networks better recognition than just black-boxes.

Classical methods for solving PDEs have been extensively developed and refined over the years, providing a basis for understanding and analyzing various physical phenomena. These methods involve discretization the PDEs using techniques as the finite difference method LeVeque (), finite element method Bathe (), finite volume method Eymard et al. () or spectral methods Trefethen (), followed by numerical solution of the resulting algebraic equations. While these methods have been successful in solving a wide range of PDEs, they often face the cure of dimensionality when parametric PDEs need to be solved in connection with optimization, optimal control, parameter identification, uncertainty quantification. The reduction of complexity for such classes of problems can be addressed with surrogate models using machine learning.

The main approaches in scientific machine learning are (i) using governing equations as loss functions with physics-informed neural networks Karniadakis et al. (); Cai et al. (); Eivazi et al. (); Raissi et al. (); (ii) learning mappings between infinite-dimensional function spaces with neural operators Li et al. (); Fanaskov and Oseledets (); Lu et al. (); Li et al. (); Tran et al. (); (iii) hybrid approaches where machine learning techniques are incorporated into classical simulations Brunton and Kutz (); Schnell and Thuerey (); Hsieh et al. (); Ingraham et al. ().

These surrogate models have shown significant potential in solving parametric PDEs, but a critical aspect of their development remains the availability of comprehensive datasets for validation. The accuracy and reliability of these machine learning-based approaches are highly dependent on the quality and diversity of the data used to train and test them. Without such datasets, the performance and generalization ability of these models cannot be adequately assessed, and their applicability to real-world problems may be limited. As new techniques and methods emerge in the future, the need for robust and extensive datasets will only increase. It is therefore essential to develop approaches to the curation of high quality datasets that can support the development and validation of innovative approaches to solving complex problems in different scientific and engineering domains.

With ConDiff (short for Contrast Diffusion), we propose a novel dataset that complements existing datasets and provides an approach to generating complex coefficients for parametric PDEs that can address real-world problems. We also provide a measurable metric for the complexity of the dataset and validate our approach on a list of well-known models.

The ConDiff dataset is available on the Hugging Face Hub: [https://huggingface.co/datasets/condiff/ConDiff](https://huggingface.co/datasets/condiff/ConDiff). The code with ConDiff generation, usage, validation and requirements is available at: [https://github.com/condiff-dataset/ConDiff](https://github.com/condiff-dataset/ConDiff).

## 2 ConDiff

#### Motivation

Creating a comprehensive benchmark for classes of parametric PDEs is a particular challenge for the scientific machine learning community. The main challenges in creating a comprehensive dataset are: (i) computational complexity; (ii) storage complexity for the desired dimensions of the discretized PDE and parameter space; (iii) properties of the coefficients and solution functions; (iv) relation to real-world problems. The first and second reasons illustrate a technical bottleneck in the creation of the dataset and are mostly dependent on the hardware and efficiency of the numerical method used. Properties such as coefficient smoothness, discontinuity, spatial variation of the coefficients, variance of the parametric space significantly affect the complexity of the dataset and should be carefully chosen. The solution to PDEs (i.e. the ground truth for the dataset) depends on a number of numerical aspects such as choice of mesh, discretization, numerical algorithm, boundary and initial conditions. Therefore, it is very important to consider every little detail regarding different numerical schemes, PDEs, boundary and initial conditions.

Existing benchmarks and datasets cover different aspects of scientific machine learning for different classes of PDEs and can be divided into several groups. PDEBench Takamoto et al. (), PINNacle Hao et al. (), CFDBench Luo et al. () have a large number of PDEs with different boundary and initial conditions and different dimensionality and resolution. The best covered area is weather forecasting: SuperBench Ren et al. (), ClimSim Yu et al. (), DynaBench Dulny et al. (), OceanBench Johnson et al. (), ChaosBench Nathaniel et al. (). There are also domain specific datasets with applications to Lagrangian mechanics LagrangeBench Toshev et al. () and phase change phenomena BubbleML Hassan et al. (). Worth noting frameworks for differential simulations and general environments for PDEs in scientific machine learning: PDE Control Gym Bhan et al. (), PDEArena Gupta and Brandstetter (), DiffTaichi Hu et al. (), DeepXDE Lu et al. () and $\Phi_{\text{Flow}}$ Holl et al. ().

Figure: Figure 1: Visualization of the GRF (top row), the coefficient $k(x)$ generated from this GRF (middle row) and the corresponding solution of the equation (1) (bottom row) for a sampled PDEs with grid $128\times 128$ and $\sigma^{2}=2.0$.
Refer to caption: extracted/5650811/results/grf-coeff-sol-ind32.png

Figure: Figure 2: Cross section of the $x-$permeability field along the $z$ axis over the SPE10 model 2 with $z=4$.
Refer to caption: extracted/5650811/results/spe10_z4.png

While all of these datasets contribute significantly to the community, to the best of the authors’ knowledge there is no dataset dedicated to the very important class of academic and real-world problems, the class of parametric PDEs with random coefficients. Typically, when a new model is proposed, authors test it with a set of equations with smooth coefficients Brandstetter et al. (); Nguyen et al. (); Ripken et al. (); Bryutkin et al. (). Such coefficients do not allow important classes of industrial applications to be addressed. In section 3 we show that increasing the heterogeneity and contrast in the coefficient function leads to increasing challenges in building accurate surrogate models.

#### Problem definition

Existing benchmarks Takamoto et al. (); Hao et al. (); Luo et al. () cover a set of PDEs, both steady-state and time-dependent, with different resolutions and time lengths. In our work, we approach the problem from the other side tacking a fixed parametric PDE and generating a comprehensive set of random coefficients for it. We consider a 2D steady-state diffusion equation:

$$ $\begin{split}-&\nabla\cdot\big{(}k(x)\nabla u(x)\big{)}=f(x),~{}\text{in}~{} \Omega\\ &u(x)\Big{|}_{x\in\partial{\Omega}}=0\end{split}\,\,\,.$ (1) $$

Note that the equation (1) models not only diffusion, but also steady-state Darcy flow in porous media, steady-state heat conduction, etc. To address certain real-world problems, we use the Gaussian Random Field (GRF) to generate the field $\phi(x)$ with the following covariance models as functions of distance $d$:

- •
Cubic:
$\text{Cov}(d)=\begin{cases}\sigma^{2}\Big{(}1-7\big{(}\frac{d}{l}\big{)}^{2}+
\frac{35}{4}\big{(}\frac{d}{l}\big{)}^{3}-\frac{7}{2}\big{(}\frac{d}{l}\big{)}
^{5}+\frac{3}{4}\big{(}\frac{d}{l}\big{)}^{7}\Big{)}\,,&d<l\\
0\,,&d\geq l\end{cases}\,.$
(2)
- •
Exponential:
$\text{Cov}(d)=\sigma^{2}\exp{\Big{(}-\frac{d}{l}\Big{)}}\,.$
(3)
- •
Gaussian:
$\text{Cov}(d)=\sigma^{2}\exp{\Big{(}-\frac{d^{2}}{l^{2}}\Big{)}}\,.$
(4)

The correlation length in each dataset is $l=0.05$ and the complexity of a resulting dataset is controlled by variance $\sigma^{2}$. The forcing term $f(x)$ is sampled from the standard normal distribution for each sampled PDE in each dataset. The resulting coefficient $k(x)$ is obtained with:

$$ $k(x)=\exp{\big{(}\phi(x)\big{)}}\,.$ (5) $$

We propose to measure the complexity of the generated GRF with the global contrast in the field $\phi(x)$:

$$ $\text{contrast}=\exp\Big{(}\max\big{(}\phi(x)\big{)}-\min\big{(}\phi(x)\big{)} \Big{)}\,.$ (6) $$

#### Complexity grows with variance

By increasing the variance $\sigma^{2}$ one can obtain a higher contrast (6) and thus a higher complexity of the PDE. This is a well-known phenomenon in applied numerical analysis and can be easily observed empirically. We illustrate this behaviour with the condition number $\kappa(A)$ of the matrices $A$ obtained with discretization of the equation (1).

In the Table 1 one can observe that increasing $\sigma^{2}$ leads to a higher condition number $\kappa(A)=|\lambda_{\max}|\big{/}|\lambda_{\min}|$ of the discretized differential operator Capizzano (). The condition number is closely related to the performance of the numerical methods used to solve PDEs Benzi et al. (); Elman et al. (). A high condition number indicates that small changes in the input can lead to large changes in the output, making the problem ill-conditioned. This is particularly important in PDEs, where small perturbations can significantly affect the solution. Also, if iterative methods are used to solve the discretized PDE, a larger condition number means a larger number of iterations for unpreconditioned and most of preconditioned iterative methods Saad ().

#### Connection to real-world problems

All of the above reasoning is done with regard to the frequent occurrence of such tasks in real world Hashmi (); Massimo (); Carr and Turner (); Oristaglio and Hohmann (); Muravleva et al. (), including composite materials modeling, heat transfer, geophysical problems, fluid flow modeling. In Figure 2 one can see a cross section of the $x-$permeability field along the $z$ axis over the SPE10 model 2 benchmark Christie and Blunt (). The term permeability is used to denote the coefficients of the above equation when considering flow in porous media. This field is very similar to the ConDiff samples in Figure 1.

This benchmark is well known in the field of reservoir modelling and fluid flow in porous media. SPE10 model 2 poses a significant challenge for the tasks of uncertainty quantification, upscaling and multiphase fluid flow modelling.

**Table 1: Summary of the ConDiff with min, mean and max values of the contrast (6). ^1Condition number $\kappa(A)$ is calculated for a single sampled discretized (1).**
| Covariance | Variance | Min contrast | Mean contrast | Max contrast | $\kappa^{1}(A)$ |
| --- | --- | --- | --- | --- | --- |
| $\text{Grid}~{}64\times 64$ |  |  |  |  |  |
| Cubic | $0.1$ | $7.0\cdot 10^{0}$ | $1.0\cdot 10^{1}$ | $1.5\cdot 10^{1}$ | $3.6\cdot 10^{3}$ |
| $0.4$ | $5.0\cdot 10^{1}$ | $9.6\cdot 10^{1}$ | $2.5\cdot 10^{2}$ | $7.3\cdot 10^{3}$ |  |
| $1.0$ | $6.0\cdot 10^{2}$ | $8.3\cdot 10^{2}$ | $1.0\cdot 10^{3}$ | $2.0\cdot 10^{4}$ |  |
| $2.0$ | $8.0\cdot 10^{4}$ | $8.9\cdot 10^{4}$ | $1.0\cdot 10^{5}$ | $1.8\cdot 10^{5}$ |  |
| Exp | $0.1$ | $6.0\cdot 10^{0}$ | $9.0\cdot 10^{0}$ | $1.5\cdot 10^{1}$ | $4.3\cdot 10^{3}$ |
| $0.4$ | $5.0\cdot 10^{1}$ | $8.5\cdot 10^{1}$ | $2.3\cdot 10^{2}$ | $5.2\cdot 10^{3}$ |  |
| $1.0$ | $6.0\cdot 10^{2}$ | $7.9\cdot 10^{2}$ | $1.0\cdot 10^{3}$ | $1.7\cdot 10^{4}$ |  |
| $2.0$ | $8.0\cdot 10^{4}$ | $8.9\cdot 10^{4}$ | $1.0\cdot 10^{5}$ | $1.9\cdot 10^{5}$ |  |
| Gauss | $0.1$ | $5.0\cdot 10^{0}$ | $8.0\cdot 10^{0}$ | $1.4\cdot 10^{1}$ | $4.1\cdot 10^{3}$ |
| 0$.4$ | $5.0\cdot 10^{1}$ | $7.5\cdot 10^{1}$ | $2.3\cdot 10^{2}$ | $8.1\cdot 10^{3}$ |  |
| $1.0$ | $6.0\cdot 10^{2}$ | $7.7\cdot 10^{2}$ | $1.0\cdot 10^{3}$ | $2.4\cdot 10^{4}$ |  |
| $2.0$ | $8.0\cdot 10^{4}$ | $8.9\cdot 10^{4}$ | $1.0\cdot 10^{5}$ | $8.8\cdot 10^{5}$ |  |
| $\text{Grid}~{}128\times 128$ |  |  |  |  |  |
| Cubic | $0.1$ | $8.0\cdot 10^{0}$ | $1.1\cdot 10^{1}$ | $1.5\cdot 10^{1}$ | $1.6\cdot 10^{4}$ |
| $0.4$ | $5.5\cdot 10^{1}$ | $1.3\cdot 10^{2}$ | $2.5\cdot 10^{2}$ | $3.8\cdot 10^{4}$ |  |
| $1.0$ | $6.0\cdot 10^{2}$ | $8.8\cdot 10^{2}$ | $1.0\cdot 10^{3}$ | $1.0\cdot 10^{5}$ |  |
| $2.0$ | $8.0\cdot 10^{4}$ | $8.9\cdot 10^{4}$ | $1.0\cdot 10^{5}$ | $1.2\cdot 10^{6}$ |  |
| Exp | $0.1$ | $6.0\cdot 10^{0}$ | $1.0\cdot 10^{1}$ | $1.5\cdot 10^{1}$ | $1.7\cdot 10^{4}$ |
| $0.4$ | $5.1\cdot 10^{1}$ | $1.1\cdot 10^{2}$ | $2.5\cdot 10^{2}$ | $3.3\cdot 10^{4}$ |  |
| $1.0$ | $6.0\cdot 10^{2}$ | $8.3\cdot 10^{2}$ | $1.0\cdot 10^{3}$ | $9.7\cdot 10^{4}$ |  |
| $2.0$ | $8.0\cdot 10^{4}$ | $8.9\cdot 10^{4}$ | $1.0\cdot 10^{5}$ | $6.3\cdot 10^{5}$ |  |
| Gauss | $0.1$ | $5.0\cdot 10^{0}$ | $8.0\cdot 10^{0}$ | $1.4\cdot 10^{1}$ | $1.8\cdot 10^{4}$ |
| $0.4$ | $5.0\cdot 10^{1}$ | $7.8\cdot 10^{1}$ | $2.5\cdot 10^{2}$ | $7.2\cdot 10^{4}$ |  |
| $1.0$ | $6.0\cdot 10^{2}$ | $7.7\cdot 10^{2}$ | $1.0\cdot 10^{3}$ | $1.6\cdot 10^{5}$ |  |
| $2.0$ | $8.0\cdot 10^{4}$ | $8.9\cdot 10^{4}$ | $1.0\cdot 10^{5}$ | $1.5\cdot 10^{6}$ |  |

#### Dataset description

To generate the fields $\phi(x)$ we use the highly efficient parafields library(^1^11https://github.com/parafields/parafields) with C++ backend. We use covariance models from $\{\text{cubic},\text{exponential},\text{Gaussian}\}$ with 4 variance values from $\{0.1,0.4,1.0,2.0\}$. We use the forcing term $f(x)\sim\mathcal{N}(0,1)$. A Dirichlet boundary condition is set for each coefficient realization since boundary conditions do not contribute significantly to the resulting complexity Capizzano (). The ground truth solution is obtained using cell-centered second-order finite volume method. The coefficients are in the center of cells, the values are in the nodes.

For each parameter set, we generate $1000$ training and $200$ test realizations of the diffusion equation (1) on $64\times 64$ and $128\times 128$ grids. Note that datasets with the same field parameters but different grid sizes are generated independently and do not represent the same field. The fixed geometry of ConDiff allows PDEs with different fields $\phi(x)$ to be compared without fear that different geometries will interfere with a fair comparison across different coefficient functions. To control the complexity of the generated PDEs realizations, we set contrast bounds during generation as follows:

- •
$\sigma^{2}=0.1,~{}\text{contrast}\in[5,15]$,
- •
$\sigma^{2}=0.4,~{}\text{contrast}\in[50,250]$,
- •
$\sigma^{2}=1.0,~{}\text{contrast}\in[6\cdot 10^{2},10^{3}]$,
- •
$\sigma^{2}=2.0,~{}\text{contrast}\in[8\cdot 10^{4},10^{5}]$.

In total, ConDiff consists of 24 PDEs with different GRFs and grid sizes. Table 1 summarizes the properties of ConDiff. Figure 3 illustrates the contrast distributions. Coming back to the permeability cross section of SPE10 model 2 (Figure 2), it has $\text{contrast}=2.5\cdot 10^{6}$ according to (6). We want to emphasize that although the most complex coefficient of ConDiff is smaller by an order of magnitude compared to the cross section of SPE10 model 2, our experiments show that this coefficient is too complex for the chosen models to predict well.

Figure: Figure 3: GRF contrast distribution for PDEs from Table 1.
Refer to caption: extracted/5650811/results/distributions_datasets.png

## 3 Experiments

#### Models

We do not attempt to benchmark every scientific machine learning surrogate model on the ConDiff. We have selected the following list of models to validate on the ConDiff: Spectral Neural Operator (SNO) Fanaskov and Oseledets (), Factorized Fourier Neural Operator (F-FNO) Tran et al. (), Dilated ResNet (DilResNet) Yu et al. () and U-Net Ronneberger et al. (). Neural operators FNO and SNO are both types of neural networks designed to learn mappings between function spaces, in particular to solve PDEs. Neural operators are designed to be universal approximators of continuous operators acting between Banach spaces and to be discretization invariant, meaning that they can handle different discretizations of the underlying function spaces without requiring changes to the model. DilResNet and U-Net are classical neural network models originating from the field of computer vision (CV). Both models have shown their applicability beyond CV and have been used extensively for modeling physical phenomena Stachenfeld et al. (); Ma et al. (). More details about the models used can be found in the Appendix A.

#### Experiment environment

For training neural networks we use frameworks from the JAX Bradbury et al. () ecosystem: Equinox Kidger and Garcia () and Optax DeepMind et al. (). The loss function used is the relative $L_{2}$ loss:

$$ $L_{2}=\frac{1}{N}\sum_{i=1}^{N}\frac{\|\hat{y_{i}}-y_{i}\|_{2}}{\|y_{i}\|_{2}} ~{}.$ (7) $$

Training samples for the models are the values of the coefficient function $k(x)$ and the forcing term $f(x)$ in the grid cells. Targets are the values of the solution function $u(x)$ in the grid cells. We also use (7) as a primary performance metric, assessing the quality of the models’ predictions, and report averaged values over the test set with standard deviation.

For all the problems we train for $400$ epochs for $\text{grid}=64$ and for $500$ epochs for $\text{grid}=128$. We use the AdamW optimizer with an initial learning rate equals to $10^{-3}$ and a weight decay equals to $10^{-2}$. We use a learning rate schedule that halves the learning rate every $50$ epochs. Each PDE realization has a dataset size of $1000$ training samples and $200$ test samples. We use a single GPU Nvidia Tesla V100 16Gb for training on $\text{grid}=64$ and a single GPU Nvidia A40 48Gb for training on $\text{grid}=128$.

**Table 2: Results for Poisson equation.**
| Grid | SNO | F-FNO | DilResNet | U-Net |
| --- | --- | --- | --- | --- |
| 64 | $0.056\pm 0.018$ | $0.027\pm 0.008$ | $0.018\pm 0.005$ | $0.020\pm 0.007$ |
| 128 | $0.073\pm 0.021$ | $0.047\pm 0.013$ | $0.063\pm 0.016$ | $0.267\pm 0.049$ |

**Table 3: Performance comparison of the models on the PDEs with the $64\times 64$ grid from ConDiff.**
| Covariance | Variance | SNO | F-FNO | DilResNet | U-Net |
| --- | --- | --- | --- | --- | --- |
| Cubic | $0.1$ | $0.09\pm 0.02$ | $0.07\pm 0.02$ | $0.07\pm 0.02$ | $0.08\pm 0.02$ |
| $0.4$ | $0.15\pm 0.04$ | $0.14\pm 0.03$ | $0.14\pm 0.03$ | $0.17\pm 0.04$ |  |
| $1.0$ | $0.23\pm 0.06$ | $0.22\pm 0.06$ | $0.22\pm 0.06$ | $0.24\pm 0.06$ |  |
| $2.0$ | $0.35\pm 0.10$ | $0.34\pm 0.09$ | $0.35\pm 0.10$ | $0.42\pm 0.10$ |  |
| Exp | $0.1$ | $0.12\pm 0.03$ | $0.11\pm 0.03$ | $0.11\pm 0.03$ | $0.12\pm 0.04$ |
| $0.4$ | $0.21\pm 0.06$ | $0.21\pm 0.06$ | $0.20\pm 0.06$ | $0.26\pm 0.07$ |  |
| $1.0$ | $0.33\pm 0.09$ | $0.34\pm 0.09$ | $0.36\pm 0.09$ | $0.35\pm 0.09$ |  |
| $2.0$ | $0.59\pm 0.14$ | $0.58\pm 0.14$ | $0.60\pm 0.13$ | $0.64\pm 0.13$ |  |
| Gauss | $0.1$ | $0.12\pm 0.04$ | $0.11\pm 0.04$ | $0.11\pm 0.03$ | $0.12\pm 0.03$ |
| $0.4$ | $0.23\pm 0.06$ | $0.22\pm 0.06$ | $0.21\pm 0.06$ | $0.25\pm 0.06$ |  |
| $1.0$ | $0.38\pm 0.08$ | $0.37\pm 0.09$ | $0.38\pm 0.09$ | $0.39\pm 0.09$ |  |
| $2.0$ | $0.66\pm 0.14$ | $0.65\pm 0.14$ | $0.66\pm 0.13$ | $0.72\pm 0.24$ |  |

**Table 4: Performance comparison of SNO and F-FNO on the PDEs with the $128\times 128$ grid from ConDiff.**
| Covariance | Variance | SNO | F-FNO |
| --- | --- | --- | --- |
| Cubic | $0.1$ | $0.09\pm 0.03$ | $0.08\pm 0.02$ |
| $0.4$ | $0.15\pm 0.04$ | $0.14\pm 0.04$ |  |
| $1.0$ | $0.23\pm 0.06$ | $0.22\pm 0.06$ |  |
| $2.0$ | $0.36\pm 0.11$ | $0.36\pm 0.10$ |  |
| Exp | $0.1$ | $0.13\pm 0.03$ | $0.12\pm 0.03$ |
| $0.4$ | $0.21\pm 0.07$ | $0.21\pm 0.06$ |  |
| $1.0$ | $0.33\pm 0.09$ | $0.33\pm 0.08$ |  |
| $2.0$ | $0.58\pm 0.15$ | $0.57\pm 0.13$ |  |
| Gauss | $0.1$ | $0.13\pm 0.04$ | $0.12\pm 0.03$ |
| $0.4$ | $0.23\pm 0.06$ | $0.23\pm 0.06$ |  |
| $1.0$ | $0.37\pm 0.10$ | $0.37\pm 0.10$ |  |
| $2.0$ | $0.68\pm 0.13$ | $0.66\pm 0.13$ |  |

#### Validation on ConDiff

We start the experiments with the Poisson equation and consider it as a special case of (1) with $k(x)=1$ and $\text{contrast}=1$. All models achieve an accuracy of the order of $10^{-2}$ (Table 2). Increasing the grid size leads to moderate increases in error, except for the U-Net for which the error increases by an order of magnitude.

The diffusion equation for grid $64$ (Table 3) with covariances (2), (3) and (4) are more challenging for the models. While the performance on the diffusion equation with cubic covariance with $\sigma^{2}=0.1$ is comparable to the performance on the Poisson equation, the error on the diffusion equation with exponential and Gaussian covariances is already an order of magnitude higher. Increasing $\sigma^{2}$ leads to worse performance of each model on each PDE. The most complex PDE is the one generated with the Gaussian covariance model in GRF, which is also consistent with the condition number estimation in Table 1. Interestingly, the performance of FNO and SNO models on PDEs with grid $128$ is not much different from PDEs on grid $64$ (Table 4).

**Table 5: Generalization of the models to unseen PDEs with different GRF covariance model with $64\times 64$ grid and $\sigma^{2}=0.1$.**
|  | SNO | F-FNO |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $\text{Train}\big{\backslash}\text{Test}$ | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.09\pm 0.02$ | $0.12\pm 0.04$ | $0.12\pm 0.03$ | $0.07\pm 0.02$ | $0.11\pm 0.03$ | $0.11\pm 0.03$ |
| Exp | $0.09\pm 0.03$ | $0.12\pm 0.03$ | $0.12\pm 0.04$ | $0.08\pm 0.03$ | $0.11\pm 0.03$ | $0.11\pm 0.04$ |
| Gauss | $0.09\pm 0.03$ | $0.12\pm 0.03$ | $0.12\pm 0.03$ | $0.08\pm 0.02$ | $0.11\pm 0.03$ | $0.11\pm 0.04$ |
|  | DilResNet | U-Net |  |  |  |  |
|  | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.07\pm 0.02$ | $0.11\pm 0.04$ | $0.11\pm 0.03$ | $0.08\pm 0.02$ | $0.12\pm 0.03$ | $0.12\pm 0.03$ |
| Exp | $0.07\pm 0.02$ | $0.11\pm 0.03$ | $0.11\pm 0.03$ | $0.08\pm 0.03$ | $0.11\pm 0.03$ | $0.11\pm 0.04$ |
| Gauss | $0.17\pm 0.06$ | $0.25\pm 0.09$ | $0.11\pm 0.04$ | $0.08\pm 0.02$ | $0.12\pm 0.04$ | $0.12\pm 0.03$ |

**Table 6: Generalization of the models to unseen PDEs with different GRF covariance model with $64\times 64$ grid and $\sigma^{2}=0.4$.**
|  | SNO | F-FNO |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $\text{Train}\big{\backslash}\text{Test}$ | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.15\pm 0.04$ | $0.22\pm 0.06$ | $0.22\pm 0.07$ | $0.14\pm 0.03$ | $0.21\pm 0.06$ | $0.21\pm 0.07$ |
| Exp | $0.18\pm 0.05$ | $0.21\pm 0.06$ | $0.22\pm 0.06$ | $0.15\pm 0.04$ | $0.21\pm 0.06$ | $0.22\pm 0.07$ |
| Gauss | $0.17\pm 0.05$ | $0.22\pm 0.06$ | $0.23\pm 0.07$ | $0.15\pm 0.04$ | $0.21\pm 0.07$ | $0.22\pm 0.06$ |
|  | DilResNet | U-Net |  |  |  |  |
|  | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.14\pm 0.04$ | $0.23\pm 0.07$ | $0.23\pm 0.07$ | $0.17\pm 0.06$ | $0.24\pm 0.07$ | $0.24\pm 0.07$ |
| Exp | $0.14\pm 0.04$ | $0.20\pm 0.06$ | $0.22\pm 0.06$ | $0.23\pm 0.08$ | $0.26\pm 0.07$ | $0.27\pm 0.08$ |
| Gauss | $0.30\pm 0.10$ | $0.24\pm 0.07$ | $0.21\pm 0.06$ | $0.21\pm 0.06$ | $0.27\pm 0.08$ | $0.26\pm 0.07$ |

**Table 7: Generalization of the models to unseen PDEs with different GRF covariance model with $64\times 64$ grid and $\sigma^{2}=1.0$.**
|  | SNO | F-FNO |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $\text{Train}\big{\backslash}\text{Test}$ | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.23\pm 0.06$ | $0.35\pm 0.09$ | $0.39\pm 0.09$ | $0.22\pm 0.06$ | $0.34\pm 0.09$ | $0.37\pm 0.09$ |
| Exp | $0.25\pm 0.06$ | $0.33\pm 0.09$ | $0.38\pm 0.09$ | $0.24\pm 0.06$ | $0.34\pm 0.09$ | $0.38\pm 0.09$ |
| Gauss | $0.24\pm 0.07$ | $0.35\pm 0.09$ | $0.38\pm 0.08$ | $0.24\pm 0.06$ | $0.35\pm 0.09$ | $0.37\pm 0.09$ |
|  | DilResNet | U-Net |  |  |  |  |
|  | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.22\pm 0.06$ | $0.35\pm 0.09$ | $0.38\pm 0.09$ | $0.24\pm 0.06$ | $0.36\pm 0.09$ | $0.38\pm 0.08$ |
| Exp | $0.25\pm 0.07$ | $0.36\pm 0.09$ | $0.38\pm 0.10$ | $0.25\pm 0.07$ | $0.35\pm 0.09$ | $0.38\pm 0.10$ |
| Gauss | $0.57\pm 0.22$ | $0.59\pm 0.22$ | $0.38\pm 0.09$ | $0.27\pm 0.07$ | $0.36\pm 0.11$ | $0.39\pm 0.09$ |

**Table 8: Generalization of the models to unseen PDEs with different GRF covariance model with $64\times 64$ grid and $\sigma^{2}=2.0$.**
|  | SNO | F-FNO |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $\text{Train}\big{\backslash}\text{Test}$ | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.35\pm 0.10$ | $0.60\pm 0.14$ | $0.70\pm 0.26$ | $0.34\pm 0.09$ | $0.61\pm 0.14$ | $0.67\pm 0.19$ |
| Exp | $0.39\pm 0.11$ | $0.59\pm 0.14$ | $0.69\pm 0.24$ | $0.39\pm 0.11$ | $0.58\pm 0.14$ | $0.66\pm 0.15$ |
| Gauss | $0.40\pm 0.11$ | $0.60\pm 0.13$ | $0.66\pm 0.14$ | $0.37\pm 0.11$ | $0.60\pm 0.13$ | $0.65\pm 0.14$ |
|  | DilResNet | U-Net |  |  |  |  |
|  | Cubic | Exp | Gauss | Cubic | Exp | Gauss |
| Cubic | $0.35\pm 0.10$ | $0.61\pm 0.14$ | $0.66\pm 0.15$ | $0.42\pm 0.10$ | $0.65\pm 0.14$ | $0.68\pm 0.14$ |
| Exp | $0.41\pm 0.10$ | $0.60\pm 0.13$ | $0.66\pm 0.17$ | $0.53\pm 0.18$ | $0.64\pm 0.13$ | $0.72\pm 0.16$ |
| Gauss | $0.72\pm 0.50$ | $0.68\pm 0.20$ | $0.66\pm 0.13$ | $0.66\pm 0.40$ | $0.69\pm 0.16$ | $0.72\pm 0.24$ |

#### Transfer between parametric spaces

Ideally, the surrogate model should handle transfers between different underlying parametric spaces of PDEs without loss of quality. In Tables 5, 6, 7, 8 show that in most experiments the error increases when training on cubic GRF and inferencing on exponential and Gaussian GRF. Conversely, the error decreases when training on Gaussian GRF and inferencing on cubic GRF.

## 4 Conclusion and Limitations

In summary, we propose a novel open-source dataset for the field of neural solving of PDEs. The unique feature of the dataset is discontinuous coefficients with high contrast for parametric PDEs from different distributions. By designing the coefficients in this way, we achieve a high complexity of the generated PDEs, which also illustrates real-world problems. The proposed complexity function allows to distinguish between the generated PDEs. We also provide code to generate new data based on the approach used in this paper. Furthermore, we validate a number of surrogate models on the ConDiff to illustrate its usefulness in the field of scientific machine learning.

It should be noted that the problems considered in this paper belong to the class of stochastic PDEs. The equation (1) has to be solved for a very large number of sampled coefficients when Monte Carlo or other methods are used to solve the stochastic PDEs. The surrogate models can help to significantly reduce the computational burden, so embedding the surrogate models tested on ConDiff into a Monte Carlo or similar stochastic PDEs solver is a reasonable next step.

Limitations of the proposed dataset are:

- 1.
For practical numerical analysis, ConDiff is generated with small and moderate variances. The case of large variances has to be studied separately.
- 2.
A linear elliptic parametric PDE is the basis of ConDiff, so other high contrast datasets are needed to test surrogate models for hyperbolic PDEs, nonlinear problems, etc.
- 3.
ConDiff is generated on a regular rectangular grid. Other meshes and geometries may be required as an evolution of ConDiff. This may require more complex computational methods to obtain the ground truth solution.
- 4.
Limited number of models validated on ConDiff. Comprehensive benchmark with different approaches (e.g. physics-informed neural networks, hybrid solvers) is needed.
- 5.
The forcing term $f(x)$ is sampled from the standard normal distributions. While in this paper we focus on the complexity arising from discontinuous coefficients with high contrast, the right-hand side of a PDE can also significantly affect the complexity of the solving PDE. The case of complex forcing terms has to be studied separately.

## Appendix A Architectures

In this section, we discuss the architectures used in more detail and provide information on the training procedures and hyperparameters used. The list of used models is:

- 1.
F-FNO – Factorized Fourier Neural Operator (F-FNO) from Tran et al. [].
- 2.
fSNO – Spectral Neural Operator (SNO). The construction mirrors FNO, but instead of FFT, a transformation based on Gauss quadratures is used Fanaskov and Oseledets [].
- 3.
DilResNet – Dilated Residual Network from Yu et al. [], Stachenfeld et al. [].
- 4.
U-Net – classical computer vision architecture introduced in Ronneberger et al. [].

### A.1 F-FNO

Unlike the original Li et al. [], the authors of Tran et al. [] proposed to changing the operator layer to:

$$ $z^{\ell+1}=z^{\ell}+\sigma\Big{[}W_{2}^{(\ell)}\sigma\Big{(}W_{1}^{(\ell)} \mathcal{K}^{(\ell)}(z^{(\ell)})+b_{1}^{(\ell)}\Big{)}+b_{2}^{(\ell)}\Big{]},$ $$

where $\sigma$ is an activation function, $W_{1}$ and $W_{2}$ are weight matrices in the physical space, $b_{1}$ and $b_{2}$ are bias vectors and

$$ $\mathcal{K}^{(\ell)}\big{(}z^{(\ell)}\big{)}=\sum_{d\in D}\Big{[}\text{IFFT} \big{(}R_{d}^{(\ell)}\cdot\text{FFT}_{d}(z^{\ell})\big{)}\Big{]},$ $$

where $R_{d}$ is a Fourier domain weight matrix, FFT and IFFT are Fast Fourier and inverse Fast Fourier transforms.

F-FNO has an encoder-processor-decoder architecture. We used the following parameters: $4$ Fourier layers in the processor, $12$ modes and GeLU as the activation function. We used $48$ features in the processor.

### A.2 SNO

We utilized spectral neural operators (SNO) Fanaskov and Oseledets [] with linear integral kernels:

$$ $u\leftarrow\int dxA_{ij}p_{j}(x)\left(p_{i},u\right)~{},$ $$

where $p_{j}(x)$ are orthogonal or trigonometric polynomials.

These linear integral kernels are an extension of the integral kernels used in the FNO Li et al. []. More specifically, starting from the input function $u^{n}$, we produce the output function $u^{n+1}$, which is later transformed by nonlinear activation. The transformation depends on the set of polynomials $p_{j}$ that form a suitable basis for the problem at hand (e.g. trigonometric polynomials, Chebyshev polynomials, etc.). These polynomials are chosen beforehand and do not change during training. The transformation is naturally divided into three parts: analysis, processing, synthesis.

At the analysis stage, we find a discrete representation of the input function by projecting it onto a set of polynomials. To do this, we compute scalar products:

$$ $\alpha_{j}=\left(p_{j},u^{n}\right)=\int dxp_{j}(x)u^{n}(x)w(x)~{},$ $$

where $w(x)$ is a non-negative weight function given by the polynomial used.

At the processing stage, we process the obtained coefficients with a linear layer:

$$ $\alpha_{i}^{\textquoteleft}=\sum_{j}A_{ij}\alpha_{j}~{}.$ $$

Finally, at the synthesis stage, we recover the continuous function as the sum of the processed coefficients:

$$ $u^{n+1}=\sum_{j}p_{j}\alpha_{j}^{\textquoteleft}.$ $$

We use SNO in Fourier basis (see Fanaskov and Oseledets []) with encoder-processor-decoder architecture. The number of SNO layers is $4$ and the number of $p_{j}(x)$ is $20$. We use GeLU as activation function.

### A.3 DilResNet

The conventional dilated residual network was first proposed in Stachenfeld et al. []. In this study, the DilResNet architecture is configured with four blocks, each consisting of a sequence of convolutions with steps of $[1,2,4,8,4,2,1]$ and a kernel size of $3$. Skip connections are also applied after each block and the GeLU activation function is used.

### A.4 U-Net

We adopt the traditional U-Net architecture proposed in Ronneberger et al. []. This U-Net configuration is characterised by a series of levels, where each level has approximately half the resolution of the previous one, and the number of features is doubled. At each level, we apply a sequence of three convolutions, followed by max pooling, and then a transposed convolution for upsampling. After upsampling, three more convolutions are applied at each level. The U-Net used in this study consists of four layers and incorporates the GeLU activation function.

## Peer Reviews

confidence: 3
contribution: 1
initial_rating: 3
presentation: 3
questions: * In Tables 2–8, while the experimental results are shown for various models across different covariance functions, is there any additional insight that can be provided beyond the raw results? What is the significance of validating models using all three covariance types—cubic, exponential, and Gaussian? Additionally, among the four models tested, can any conclusions be drawn regarding which model performs best?

* Typo in line 498: It seems that "44." is a typo.

* Different Input Scenarios in the Dataset: Were cases considered where both  k(x)  and f(x) are inputs, only k(x) is an input with f(x) fixed, and only f(x) is an input with k(x) fixed? Did the study explore the various cases that could be tested within the proposed dataset?
soundness: 3
strengths: Through extensive experiments, the paper demonstrates how the variance changes for each of the four models—SNO, F-FNO, DilResNet, and U-Net—across all train/test scenarios depending on differences in the covariance function.
summary: This paper proposes a challenging dataset as a benchmark for models that use neural networks as PDE solvers. The proposed dataset aims to be more representative of real-world data by incorporating discontinuous coefficients with high contrast.
weaknesses: First and foremost, I am curious about what differentiates and adds novelty to the dataset proposed in this paper compared to existing datasets created for scientific machine learning (sciML) in previous studies. As cited in the paper, various prior works, such as those by Takamoto et al. (2022), Luo et al. (2023), and Hao et al. (2023), have developed diverse types of datasets as PDE benchmarks. I wonder what major differences exist between this dataset and the ones developed in these studies. The paper mentions that previous datasets lack "a dataset dedicated to the very important class of academic and real-world problems, the class of parametric PDEs with random coefficients." I am interested in understanding more precisely how this differs from previous works. Furthermore, prior datasets include both time-dependent PDEs and steady-state solutions across various resolutions, while this dataset is constructed for a fixed parametric PDE. A more detailed explanation on the novelty of this approach would be helpful.

In relation to the discussion starting on line 159 of page 3, it would be helpful to clarify in more detail how the connection to real-world applications is established. The paper states that Figure 1, showing a dataset created with cubic, exponential, and Gaussian functions using GRF (Gaussian Random Fields), resembles the real-world data shown in Figure 2. However, it would be beneficial to explain more specifically in what ways this dataset offers a unique advantage over other existing datasets. Given the variety of datasets already available, understanding exactly how this proposed dataset provides a novel representation or closer alignment with real-world characteristics would strengthen the paper's argument.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: 1. Can you test them on recent machine learning models? for example PINN, or Causal PINN[1], and so on.
2. Can you use on more complex PDE system? 







[1] Sifan Wang, Shyam Sankaran, and Paris Perdikaris. Respecting causality is all you need for training physicsinformed neural networks. arXiv preprint arXiv:2203.07404, 2022.
soundness: 2
strengths: The authors propose a novel dataset for neural approaches to solving parametric PDEs, featuring discontinuous coefficients with high contrast across various distributions of parametric PDEs.
summary: This paper introduces a benchmark dataset designed to test the capability of machine learning models in solving complex partial differential equation (PDE) systems. The dataset is based on 2D diffusion equation with complex coefficients.
weaknesses: The PDE examined in this paper is relatively straightforward, and existing methods can solve it accurately. Additionally, the models tested on this benchmark are somewhat outdated; for instance, while the authors discuss Physics-Informed Neural Networks (PINNs), they do not include them in the experimental analysis. 

The forcing term \(f(x)\) is sampled from a standard normal distribution, keeping the focus on coefficient complexity. However, in many real-world scenarios, the complexity of the forcing term itself significantly affects PDE solutions. Investigating the effects of complex or structured forcing terms in the dataset would make ConDiff more representative of practical challenges and allow for a more comprehensive model evaluation.

The dataset is generated on a regular rectangular grid, which may not fully represent the geometrical complexity found in many real-world applications. Expanding ConDiff to include different grid types or irregular geometries, although computationally demanding, would improve its relevance for practical applications involving complex domains.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: - line 48-62: there seems to be some confusion about the different techniques between PINNs, FNO,... as these have very different purposes (e.g. solving the PDE by minimizing its residual, or learning the solution operator from input-output data pairs).
- line 259: I'm not sure if it's correct to say that the boundary condition does not influence the solution's regularity behavior
- line 262: how do you validate the method, what is the error of the finite volume method? One could use a method of manifactured solution
- line 322: This should be an L^2 loss and include the quadrature weights. In this case it cancels out because of the uniform discretization.
- Tables 2-4: what does that matter?
- line 480: I believe it is neither possible to reach machine-precision accuracy (due to optimization) nor desirable (to avoid overfitting). Additionally, the solution in the benchmark are not resolved to machine precision.
soundness: 2
strengths: One of the issue in the field of operator learning is the lack of availabiliy of benchmarks, making the comparison between the profusion of neural operator architectures challenging.
This paper aims to address this challenge by introducing a dataset that one can easily use to compare against other techniques.
summary: This work introduces a dataset called ConDiff for benchmarking operator learning techniques, such as the Fourier Neural Operator (FNO) and DeepONets. The authors consider a two-dimensional diffusion equation of the form $-\nabla\cdot(k(x)\nabla u) = f$ on a square domain with homogeneous Dirichlet boundary conditions and generate different datasets of the form $(k_j,u_j)_{j=1}^N$, where the diffusion coefficients $k_j$ are sampled from Gaussian processes with different covariance kernels. The chosen kernels include cubic, exponential, and Gaussian, with a range of lengthscale parameters.
weaknesses: There are several weaknesses in this work, summarized in the following points:
1. lack of physical relevance. The authors mention that the class of problems is of "great academic interest" and "shortens the gap with real-world problems". However, they only consider one single simple linear diffusion equation.
2. the authors only consider a simple square domain with rectangular grid. One would typically want to evaluate whether certain architectures can deal with non-periodic boundary condition or non-structured data, which would be closer to real-world applications. This is mentionned in point 3 of the limitations but I do not see why one cannot generate such datasets using a finite element method.
3. In the discussion, the authors state that the "unique feature of the dataset is discontinuous coefficients". However, this has been previously considered in the original FNO paper, and the dataset is publicly available online (containing also more challenging problems), hence the lack of originality of the paper.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 3
questions: "For all the problems we train for 400 epochs for grid = 64 and for 500 epochs for grid = 128."
- why not have more grid resolution, e.g. 32, 256, 512?

Limitations:
"The forcing term f(x) is sampled from the standard normal distributions. While in this pa-
per we focus on the complexity arising from discontinuous coefficients with high contrast,
the right-hand side of a PDE can also significantly affect the complexity of the solving
PDE. The case of complex forcing terms has to be studied separately."
- This severely limits the type of solutions, since the forcing term is discontinuous (noise).
Why not also use a GFF to generate smoother forcing terms?  Or use step functions with some geometry, e.g. piecewise rectangles?
soundness: 3
strengths: The paper indeed generates benchmark data for the diffusion equation.
summary: This is a benchmark and dataset paper which presents a new dataset for scientific ML.  ConDiff focuses on the diffusion equation with spatially dependent coefficients, in particular, discontinuous coefficients.   It provides a large number of problem instances, for a single equation, in contrast to other benchmarks which provide fewer instances for many equations.
weaknesses: There are already a large number of benchmark dataset for scientific machine learning.
Lines 90-11 mention about 16 of them.  It's not clear how significant another one will be. 
This benchmark focuses on just a single equation - but this equation is already in other datasets.  So the novelty is mainly in the choice of parameters of the equation.  This is a rather minor contribution.  The existing  benchmarks aim to be either: very broad, to compare performance of solvers which can handle many different equations, or very applications specific, e.g. specific and realistic physics problems.  The diffusion equation is a foundational equation, but it's is already treated in the benchmarks.  If a user want to dive deeper into this particular equation, it's quite easy to generate solutions using the same methodology as proposed in here: use a numerical solver with a parameterized set of coefficients (e.g. GFF) to generate solutions.

Meta-review: This paper introduces ConDiff, a new benchmark dataset designed to evaluate machine learning models, particularly neural network-based PDE solvers like FNOs and DeepONets, in solving the 2D diffusion equation with spatially dependent, discontinuous coefficients. It provides a large number of problem instances with varying diffusion coefficients sampled from Gaussian processes with different covariance kernels (cubic, exponential, and Gaussian) and length-scale parameters.

Reviewers questioned the physical relevance of solving such a simple PDE with a trivial square domain, for which traditional methods based on DG and multigrid techniques are already optimal and which ML techniques routinely solve efficiently. Furthermore, it is unclear what added value this dataset offers compared to other benchmarks that also contain 2D diffusion data. Consequently, the contribution of this dataset is minimal; therefore, I recommend rejection.

Decision: Reject