Metamizer: A Versatile Neural Optimizer for Fast and Accurate Physics Simulations

## Abstract

Abstract Efficient physics simulations are essential for numerous applications, ranging from realistic cloth animations or smoke effects in video games, to analyzing pollutant dispersion in environmental sciences, to calculating vehicle drag coefficients in engineering applications. Unfortunately, analytical solutions to the underlying physical equations are rarely available, and numerical solutions require high computational resources. Latest developments in the field of physics-based Deep Learning have led to promising efficiency improvements but still suffer from limited generalization capabilities and low accuracy compared to numerical solvers. In this work, we introduce Metamizer, a novel neural optimizer that iteratively solves a wide range of physical systems with high accuracy by minimizing a physics-based loss function. To this end, our approach leverages a scale-invariant architecture that enhances gradient descent updates to accelerate convergence. Since the neural network itself acts as an optimizer, training this neural optimizer falls into the category of meta-optimization approaches. We demonstrate that Metamizer achieves unprecedented accuracy for deep learning based approaches - sometimes approaching machine precision - across multiple PDEs after training on the Laplace, advection-diffusion and incompressible Navier-Stokes equation as well as on cloth simulations. Remarkably, the model also generalizes to PDEs that were not covered during training such as the Poisson, wave and Burgers equation. Our results suggest that Metamizer could have a profound impact on future numerical solvers, paving the way for fast and accurate neural physics simulations without the need for retraining.

## 1 Introduction

Countless physical systems can be described by partial differential equations (PDEs). For example, electrostatic or gravitational fields can be described by the Poisson equation, heat diffusion or pollutant dispersion can be described by the advection-diffusion equation, pressure waves by the wave equation, fluids by the incompressible Navier-Stokes equation, and so on.
Unfortunately, most of these equations do not possess analytical solutions and numerical solutions are computationally expensive.
To reduce the computational burden of numerical methods, numerous neural surrogate models have been developed in recent years that achieve fast and highly efficient physics simulations (Thuerey et al., ; Cuomo et al., ). However, the accuracy of current neural surrogates still remains limited compared to numerical solvers that can iteratively improve a solution to a high desired accuracy. In addition, neural surrogate models are usually tailored to specific PDEs and do not generalize well across multiple different PDEs without retraining. Zero-shot PDE solvers, analogous to current endeavours in foundational models for natural language processing (Brown, ) or computer vision (Radford et al., ; Kirillov et al., ), could be of great benefit for instant realisitc animations in computer games or computer aided engineering.

Thus, in this work, we propose Metamizer, a neural optimizer to accelerate gradient descent for fast and accurate physics simulations.
To this end, we introduce a novel *scale-invariant architecture* that suggests improved gradient descent update steps that are independent of arbitrary scalings of the loss function or loss domain.
We train Metamizer *without any training data* directly on the physics-based loss that it is supposed to optimize.
After training, unlike traditional gradient descent methods, Metamizer requires *no hyperparameters* such as learning rate or momentum, and allows *tradeoffs between runtime and accuracy* for *fast and highly accurate* results.
We evaluate the very same network (same architecture and same weights) across a wide range of linear and non-linear PDEs and showcase its *generalization capabilites - even to PDEs that were not covered during training*.

Code as well as a pretrained model will be made available upon acceptance.

## 2 Related Work

##### Numerical Methods

that do not rely on deep learning to solve PDEs for example by means of finite differences or finite elements require high computational resources (Farmaga et al., ).
Thus, already in the last century several specialized methods have been developed to accelerate for example fluid simulations (Chen & Doolen, ; Stam, ) or cloth simulations (Baraff & Witkin, ). Unfortunately, such specialized methods do not generalize across various different PDEs as they exploit domain-specific assumptions and are still relatively expensive to compute.

##### Gradient Descent

based methods like AdaGrad (Duchi et al., ), Adam (Kingma & Ba, ; Reddi et al., ), AdamW (Loshchilov, ) and many more (Norouzi & Ebrahimi, ) build the backbone of most Deep Learning approaches. These optimization methods can be generally applied to a wide range of differentiable minimization problems including solving PDEs (Nurbekyan et al., ). However, Gradient Descent methods require tuning of hyperparameters such as learning rate or momentum and
many iterations until convergence. As a result, they are not yet efficient enough for real-time physics simulations.

##### Meta-Learning

, and in particular ”Learning to Optimize” (L2O), is the field of research that deals with the automatic improvement of optimization algorithms, typically by means of machine learning techniques.
Pytorch libraries such as ”higher” (Grefenstette et al., ) or ”learn2learn” (Arnold et al., ) allow for example to automatically optimize hyperparameters of optimizers like ADAM through gradient descent.
Andrychowicz et al. () train LSTMs to optimize neural networks for classification tasks and style transfer
and Chen et al. () present several training techniques that suggest there is still room for potential improvement in the field L2O.
However, to the best of our knowledge, L2O has not yet been applied in the context of physics simulations.

##### Data-driven Deep Learning approaches

have been widely established to generate neural surrogate models for efficient physics simulations for example in context of fluid dynamics (Tompson et al., ; Sanchez-Gonzalez et al., ), cloth simulations (Pfaff et al., ), weather forecasting (Bonev et al., ; Lam et al., ) or aerodynamics (Li et al., ).
While these methods allow efficient simulations at coarse spatial and temporal resolutions or on physical systems without knowledge of the underlying mathematical equations, they require large amounts of high-quality training data, which can be expensive to generate.
Furthermore, generalization beyond the training data is often fairly limited.

##### Physics-informed Deep Learning approaches

as introduced by Karniadakis et al. () have become a new and quickly growing field of research (Cuomo et al., ) to solve PDEs by training an implicit neural representation with only little or no ground truth data based on a physics-informed loss. PINNs are used for a variety of problems ranging from climate simulation (Lütjens et al., ; Lai et al., ; Chen et al., ) to fluid dynamics (Cai et al., ; Ghosh et al., ), all the way to robotics (Liu et al., ; Sanyal & Roy, ; Yang et al., ).
Unfortunately, the learned implicit representations need to be retrained for every physical simulation and thus are not real-time capable.

##### Physics-constrained Deep Learning approaches

as presented by Zhu et al. () train a neural surrogate model to evolve an explicit state representation in time using a physics-constrained loss.
This method has been successfully applied to the coupled 2D Burgers equation (Geneva & Zabaras, ), incompressible fluid dynamics in 2D and 3D (Wandel et al., ; ) and cloth simulations (Santesteban et al., ; Stotko et al., ).
While these methods do not require any training data and show good generalization performance within their PDE domain, so far, they do not generalize across multiple PDEs and have limited accuracy compared to numerical solvers.

## 3 Foundations

In this chapter, we introduce some basic concepts and definitions for PDEs, explain, how finite difference schemes can be used to formulate a physics-based loss, and give some intuition to motivate our gradient-based optimization approach.

##### Partial Differential Equations (PDEs)

Partial differential equations are equations that constrain the partial derivatives of a multivariate function. In its most general form, a PDE can be written as:

$$ $F(x_{1},...,x_{n},\partial_{x_{1}}u,...,\partial_{x_{n}}u,\partial_{x_{1}}^{2} u,...)=0\,\forall x_{1},...,x_{n}\in\Omega$ (1) $$

where $u(x_{1},...,x_{n})$ is the *unknown function* that we want to solve, $x_{1},...,x_{n}\in\Omega$ are the *independent variables* inside a domain $\Omega$, and $\partial_{x_{i}}u$ denotes the partial derivative of $u$ with respect to $x_{i}$.
If one of the independent variables corresponds to time $t$, $F$ is called a *time-dependent* PDE. If $u$ does not depend on $t$, it is called *stationary*.
If $F$ is linear in $u$ and its derivatives, $F$ is called a *linear* PDE. Otherwise, $F$ is called a *non-linear* PDE.

##### Boundary and Initial Conditions

Usually, additional constrains at the domain boundaries are given such as boundary or initial conditions. Although various types of different boundary conditions could be considered, in this work, for simplicity, we focus solely on *Dirichlet boundary conditions*:

$$ $u(x_{1},...,x_{n})=d(x_{1},...,x_{n})\,\forall x_{1},...,x_{n}\in\partial\Omega$ (2) $$

Here, $d(x_{1},...,x_{n})$ directly specifies the values of $u$ at the domain boundary $\partial\Omega$.
Typically, time-dependent PDEs also require an initial state (or *initial condition*) from which the solution of a PDE evolves over time.

##### Finite Differences

To compute $u$ and its partial derivatives in a numerical manner, $u$ needs to be discretized. Here, we consider a regular grid in space (and time). This allows to compute derivatives like gradients ($\nabla$), divergence ($\nabla\cdot$), Laplace ($\Delta$) or curl ($\nabla\times$) numerically by performing convolutions with finite difference kernels.
By plugging the results of these finite difference convolutions for every grid point into the general formulation of PDEs (Equation 1), we obtain a large system of equations. If the PDE is linear, this set of equations becomes a (typically sparse) linear system of equations.

##### Time-Integration Schemes

In case of time-dependent PDEs, a strategy must be chosen to integrate the simulation in time for given initial conditions.
There exist several strategies to calculate the time derivative with finite differences:

- •
In *forward Euler* schemes, the computation of the time derivative depends only on the current time step. This explicit strategy makes computations easy and efficient, but tends to cause unstable simulation behavior for many physical systems if the time step is chosen too large. Therefore, explicit schemes are not considered in this work.
- •
In *backward Euler* schemes, the calculation of the time derivative depends on the following time step. This implicit strategy results in stable simulations, but can suffer from numerical dissipation and requires solving large systems of equations.
- •
In *Crank-Nicolson* schemes, the centered finite difference of the current and the following time step is used to calculate the time derivative. This implicit strategy results in more accurate computations compared to forward or backward Euler, but requires, like backward Euler, solving large systems of equations.

##### Physics-Based Loss

Solving large (potentially non-linear) systems of equations obtained by finite differences is a challenging task that is typically computationally expensive.
Thus, we formulate a physics-constrained loss, that penalizes the mean squared residuals (left-hand side) of Equation 1 as follows:

$$ $L(u)=\frac{1}{|\Omega|}\sum_{x_{1},...,x_{n}\in\Omega}|F(x_{1},...,x_{n}, \partial_{x_{1}}u,...,\partial_{x_{n}}u,\partial_{x_{1}}^{2}u,...)|^{2}$ (3) $$

This way, we can search for a solution $u^{*}$ of the PDE by minimizing $L$:

$$ $u^{*}=\operatorname*{arg\,min}_{u}L$ $$

Dirichlet Boundary conditions (Equation 2) can be typically applied directly by setting the corresponding boundary grid points of $u$ equal to $d$.
However, sometimes, this is not possible and we have to consider an additional loss term for the Dirichlet boundary conditions:

$$ $L_{D}(u)=\frac{1}{|\partial\Omega|}\sum_{x_{1},...,x_{n}\in\partial\Omega}|u(x _{1},...,x_{n})-d(x_{1},...,x_{n})|^{2}$ (4) $$

In this case, we have to minimize $L+L_{D}$ to solve Equation 1 and Equation 2.

##### Gradient Descent Intuition

Figure: Figure 2: Exemplary visualization of a 2D loss function with a steep valley. Naive gradient descent fails since it either converges very slowly (blue curve 1) or diverges (red curve 2). Ideally, we would like to adapt the step size and direction for accelerated convergence (green curve 3).
Refer to caption: extracted/5916378/imgs/Intuition.png

In order to minimize the physics-based loss $L$ (Equation 3 and 4), we follow a gradient descent approach.
Figure 2 provides some intuition on common problems with naive gradient descent in the vicinity of ill conditioned minima and how convergence could be improved:
If we take too small steps along the gradients (blue curve 1), convergence becomes very slow.
If we take too large steps (red curve 2), the gradient descent scheme diverges.
Ideally (green curve 3), we would first take some small initial steps to check the local gradients and avoid initial divergence. Then, the step size can be carefully increased. Instead of directly following the gradients, a better step direction can be found by taking the gradients of the previous steps into account as well. Once we come close to the minimum, the step size should be decreased again to achieve higher accuracy.
Note that scaling the loss function $L$ and thus its gradients by a constant factor $c$ does not affect its minima ($\operatorname*{arg\,min}_{u}L(u)=\operatorname*{arg\,min}_{u}c\cdot L(u)$). Thus, the gradient descent procedure should be invariant with respect to constant gradient scalings.
Similarly, scaling the coordinate system of $L$ should result in equally scaled update steps ($\operatorname*{arg\,min}_{u}L(u)=c\cdot\operatorname*{arg\,min}_{u}L(c\cdot u)$, see solid vs dashed lines in Figure 2).
This motivates our scale invariant Metamizer architecture that can deal with arbitrarily scaled gradients and can automatically adapt its step size (see Figure 3 a).

## 4 Method

Our approach aims to accelerate convergence of a gradient descent scheme by proposing improved update steps with a neural network.

##### Metamizer Architecture

Figure: Figure 3: a) Scale invariant architecture of the neural optimizer (Metamizer). b) Training Cycle.
Refer to caption: extracted/5916378/imgs/Architecture_Cycle.png

The neural network architecture as depicted in Figure 3 a) works as follows: First, the gradient $\nabla_{u}L_{i}$ of the current iteration $i$ as well as the gradient $\nabla_{u}L_{i-1}$ and update step $\Delta u_{i-1}$ of the previous iteration $i-1$ are taken as input.
Then, the gradients of the current and previous iteration are normalized by the standard deviation of the current gradient and the step size is normalized by its standard deviation as well.
These values are then concatenated and fed into a U-Net (Ronneberger et al., ) in order to predict an update scale factor for $s_{i}$ and an update direction that gets multiplied with $s_{i}$ to obtain the update step $\Delta u_{i}$ for the next iteration.
Finally, the computed update step is applied to $u$ to obtain $u_{i+1}=u_{i}+\Delta u_{i}$.
For the initial state, we set $u_{0}=\nabla_{u}L_{0}=\Delta u_{0}=0$ and set $s_{0}=0.05$ to avoid initial divergence by making only small steps at the beginning.

These normalization and scaling measures make the network invariant with respect to arbitrary scalings of the loss or domain and allow the network to make reasonable progress at all stages of the optimization:
At the beginning, when large update steps must be made in the presence of large gradients, as well as in later optimization stages, when fine adjustments must be made in the presence of small gradients.
To allow for very high accuracies, the gradient computations and update steps need to be calculated in double precision. Nevertheless, due to the input normalizations, the CNN can be evaluated using a more efficient regular float precision.
Since the outputs of the network get scaled by $s_{i}$, the training gradients get scaled as well. This would result in huge gradients for the first optimization iterations compared to basically no training signal at later optimization iterations.
Thus, we employ gradient normalization (Chen et al., ) to normalize gradients during backpropagation and provide the network with an evenly distributed training signal during all optimization stages and across all different PDEs that we train on in parallel.
Since we want to reuse the same network for different PDEs that might also contain different numbers of channels, Metamizer optimizes all channels separately (e.g. the Laplace equation requires only 1 channel while incompressible fluids require 2 channels for the velocity and pressure field and cloth dynamics require 3 channels to describe motions in $x/y/z$ directions).

##### Training

The neural optimizer was trained in a cyclic Meta-learning manner inspired by Wandel et al. ():
As visualized in Figure 3 b), we first 0-initialize a training pool of randomized PDEs (e.g. randomized boundary conditions for the Laplace / advection-diffusion / Navier-Stokes equations). Note that in contrast to data-driven methods, no pre-generated ground truth data is needed. Then, we draw a random mini-batch that contains the current gradient $\nabla_{u}L_{i}$ of the physics-based loss as well as optimizer state information (i.e. $\nabla_{u}L_{i-1},\Delta u_{i-1},s_{i-1}$) and feed it into the Metamizer architecture to compute an update step $\Delta u_{i}$. This update step, together with updated state information (i.e. $\nabla_{u}L_{i},\Delta u_{i},s_{i}$) is fed back into the training pool to update the training environments with more realistic data and to recompute the physics-based loss $L$ based on the updated values of $u$.
This recomputed loss is then used to optimize the Metamizer with gradient descent (Adam, lr=0.001). By iterating this training cycle, the training pool becomes filled with more and more realistic training data and the neural optimizer becomes better and better at minimizing physics-based losses.
Furthermore, since this way the training pool constantly generates new training data on the fly, its memory footprint is small enough such that the entire training pool can be kept in GPU memory.
After training for ca 6 hours on a NVidia Geforce RTX 4090, we obtained a model that was able to solve various PDEs at high precision and produced the results shown in the following section.

## 5 Results

In this section, we show that the very same model (same architecture and same weights) can deal with a wide variety of PDEs - linear as well as non-linear, stationary as well as time-dependent.

### 5.1 Linear Second Order PDEs

Linear second order PDEs can be grouped into 3 categories:
Elliptic PDEs (e.g. Poisson or Laplace equation), parabolic PDEs (e.g. advection-diffusion equation) and hyperbolic PDEs (e.g. wave equation).
Here, we show that our model excels in all 3 categories:

#### 5.1.1 Poisson and Laplace equation

The Poisson equation is a stationary, elliptic PDE and is important to find e.g. minimal surfaces, steady state equilibria of the Diffusion-Equation (see Section 5.1.2) or to calculate energy potentials of electrostatic or gravitational fields (see Figure 12 in Appendix):

$$ $F=\Delta u-f=0$ (5) $$

The Laplace Equation is the homogeneous ($f=0$) special case of the Poisson equation.

Figure: Figure 4: Visualization of iterative refinements by Metamizer to solve the Laplace equation. Top row: Intermediate solutions of $u$. Middle row: corresponding gradients of the physics-based loss. Bottom row: update steps computed by Metamizer.
Refer to caption: extracted/5916378/imgs/laplace_steps.png

Figure 4 depicts intermediate solutions $u_{i}$ of the Laplace equation, the gradients of the physics-based loss $\nabla_{u}L_{i}$, and update steps $\Delta u_{i}$ computed by Metamizer at iterations 0,1,2,3,20 and 50. At the beginning, $u$ is 0-initialized.
Thus, high gradients of $L$ only appear at the domain boundaries where the boundary conditions result in discontinuities of $u$. Nevertheless, Metamizer is able to anticipate better update steps that also affect the surroundings of the boundaries to achieve faster convergence. After only 3 update steps, the shape of $u$ is visually close to the final solution. After 20 iterations, the gradients of $L$ are already in the range of $10^{-7}$ and after 50 steps, Metamizer hits the limits of machine precision (see noise in $\nabla_{u}L_{50}$ in Figure 4) rendering further improvements impossible.
Figure 5 a shows how Metamizer automatically adjusts the scaling $s_{i}$: At the very beginning, $s_{0}$ starts at a relatively small value (0.05) to avoid initial divergence. Then, for the next 5 iterations, Metamizer increases $s_{i}$ to lower $L$ more efficiently. After that, $s_{i}$ is steadily decreased for finer and finer adjustments.
Remarkably, this learned behavior corresponds very well with our ”natural” gradient descent intuition provided in Section 3.
After around 50 iterations, when machine level precision is reached, no further progress is possible and $s_{i}$ remains at a constant very low level.

We compared our approach to several other traditional gradient descent based methods (SGD, Adam, AdamW, RMSprop, Adagrad, Adadelta) and state of the art iterative sparse linear system solvers (conjugate gradients, minres, gcrotmk, gmres, lgmres).
Figure 6 shows that traditional gradient based approaches are not competitive with current sparse linear system solvers regarding runtime or accuracy. Furthermore, gradient based methods require tuning of hyperparameters such as learning rate or momentum.
On the other hand, linear system solvers are limited to linear PDEs and cannot be directly used to solve non-linear PDEs.
Metamizer maintains the generality of gradient based approaches while showing slightly better performance than linear system solvers.
Furthermore, it can be easily parallelized on GPUs which results in better scaling behavior to larger grid sizes (see also Appendix B).
On a $100\times 100$ grid, Metamizer reaches a MSR of $6.7\times 10^{-33}$ for the Laplace equation after around 50 iterations in only $0.08$ seconds.
To the best of our knowledge, there exist no other deep learning based approaches that reach this level of accuracy.

Figure: Figure 5: The scale parameter $s_{i}$ (see Figure 3 a) is automatically adjusted by Metamizer to make appropriate update steps. a) stationary Laplace Equation b) time-dependent advection-diffusion equation with 10 iterations per time step.
Refer to caption: extracted/5916378/imgs/scales.png

Figure: Figure 6: Comparison of mean squared residual errors as a function of time (left) and number of iterations (right) between Metamizer and various traditional gradient based optimizers from pytorch as well as state of the art iterative sparse linear system solvers from scipy. The target is to lower the loss in as little time as possible (bottom left corner of the left plot).
Refer to caption: extracted/5916378/imgs/comparison_iterative_sparse_solvers_2.png

#### 5.1.2 Advection-Diffusion equation

The Advection-Diffusion equation is an important parabolic PDE that describes for example heat transfer or the dispersion of pollutants. For an unknown function $u$, it can be written as:

$$ $F=\partial_{t}u-D\Delta u+\nabla\cdot(\vec{v}\,u)-R=0$ (6) $$

Here, $D$ denotes the diffusivity, $\vec{v}$ denotes the velocity field for the advection term and $R$ corresponds to the source / sink term.
In contrast to the stationary Poisson equation, the advection-diffusion equation is time-dependent.
Thus, we use the implicit backward Euler scheme for a stable unfolding of the simulation in time. For every simulation time step, we perform a certain number of optimization steps to reduce $L$.
By choosing the right number of iterations per time step we can make a trade-off between speed and accuracy.
Figure 5 b) shows how Metamizer automatically adjusts the scaling $s_{i}$ when simulating the advection-diffusion equation at 10 iterations per time step. At the beginning of each time step, the neural optimizer increases its stepsize to quickly adjust $u$. Then, $s_{i}$ is decreased again for fine adjustments.

Figure 1 shows results for a constant localized source field $R$ and a rotating velocity field $\vec{v}$ (for example a laser that heats a spinning disk). More examples with different diffusivity parameters are shown in Appendix A. Table 1 shows how accuracy increases with the number of iterations per time step.

#### 5.1.3 Wave equation

The wave equation is an important hyperbolic PDE that describes wave like phenomena such as water waves, pressure waves, electro-magnetic waves and many more. It can be written as:

$$ $F=\partial_{t}^{2}u-c^{2}\partial_{x}^{2}u=0$ (7) $$

Here, $c$ corresponds to the wave propagation speed. Figure 7 shows results for different values of $c$ that were obtained by Metamizer using an implicit Crank-Nicolson scheme. Note, that the wave equation was not included during training.
Table 1 shows accuracy levels for 1,5,20 and 100 iterations per timestep. Using only 1 iteration results in divergent simulations, 20 iterations yield fairly accurate results and at 100 iterations, the simulation runs close to machine precision.

Figure: Figure 7: Metamizer is capable of simulating waves at different propagation speeds $c$ without having seen the wave equation during training.
Refer to caption: extracted/5916378/imgs/Wave_parameters.png

### 5.2 Non-Linear PDEs

Non-linear PDEs are particularly hard to solve as the underlying systems of equations become non-linear as well. In this section, we demonstrate Metamizers capabilities to deal with non-linear PDEs.

#### 5.2.1 Incompressible Navier-Stokes equation

The incompressible Navier-Stokes equation describes the dynamics of an incompressible fluid by means of a velocity field $\vec{v}$ and a pressure field $p$:

$$ $F=\rho\left(\partial_{t}\vec{v}+\left(\vec{v}\cdot\nabla\right)\vec{v}\right)+ \nabla p-\mu\Delta\vec{v}-\vec{f}_{\textrm{ext}}=0$ (8) $$

Here, $\rho$ denotes the fluids density, $\mu$ its viscosity and $\vec{f}_{\textrm{ext}}$ external forces.
The incompressibility equation $\nabla\cdot\vec{v}=0$ can be automatically fulfilled by utilizing a vector potential $a$ and setting the velocity field $\vec{v}=\nabla\times a$.
However, using a vector potential prohibits us from directly manipulating the velocity field. Thus, we need an additional boundary loss term (Equation 4) to deal with the Dirichlet boundary conditions for the velocity field at the domain boundaries.
To compute the partial derivatives of $a,\vec{v}$ and $p$ efficiently with centered finite differences, we rely on a staggered marker and cell (MAC) grid (Harlow et al., ; Holl et al., ; Wandel et al., ).

Figure 8 shows results for various Reynolds numbers. The Reynolds number $Re$ is a unit less quantity defined as:

$$ $Re=\frac{\rho||\vec{v}||D}{\mu}$ (9) $$

where $D$ corresponds to the diameter of the obstacle and $||\vec{v}||$ to the flow velocity.
The Reynolds number has a big effect on the qualitative behavior of the flow field. For example, in case of very small Reynolds numbers, the flow-field becomes stationary and the Navier-Stokes equation can be approximated by the Stokes equation. In case of very high Reynolds numbers, the flow-field becomes turbulent and can be approximated by the Euler equation for inviscid fluids.

Figure: Figure 8: Fluids of various viscosities $\mu$ and densities $\rho$ at a wide range of Reynolds numbers can be simulated by Metamizer. The obstacle diameter is $D=10$ and the fluid velocity is $||\vec{v}||=0.5$
Refer to caption: extracted/5916378/imgs/Fluid_parameters.png

**Table 1: Mean Squared Residuals for different time-dependent PDEs at different iterations per time step**
| PDE \$\frac{\textrm{Iterations}}{\textrm{Time Step}}$ | 1 | 5 | 20 | 100 | $\frac{\textrm{Iterations}}{\textrm{Second}}$ |
| --- | --- | --- | --- | --- | --- |
| Advection-Diffusion | $7.9\times 10^{-5}$ | $2.0\times 10^{-7}$ | $5.6\times 10^{-11}$ | $2.2\times 10^{-33}$ | 420 |
| Wave | $4.4\times 10^{24}$ | $4.6\times 10^{-4}$ | $3.2\times 10^{-8}$ | $3.5\times 10^{-31}$ | 420 |
| Navier-Stokes $L_{p}$ | $7.4\times 10^{-6}$ | $2.0\times 10^{-5}$ | $2.4\times 10^{-8}$ | $2.4\times 10^{-16}$ | 230 |
| $L_{d}$ | $1.7\times 10^{-7}$ | $3.7\times 10^{-7}$ | $9.5\times 10^{-9}$ | $1.9\times 10^{-17}$ |  |
| 2D coupled Burgers | $4.1\times 10^{-3}$ | $8.5\times 10^{-6}$ | $2.5\times 10^{-10}$ | $5.0\times 10^{-32}$ | 330 |

Table 1 shows results for $L_{p}$ (mean squared residuals of Equation 8) and $L_{d}$ (mean squared residuals of $\nabla\cdot\vec{v}$ that capture if the Dirichlet boundary conditions are not properly met).
$L_{p}$ and $L_{d}$ are lower for only 1 iteration compared to 5 iterations per time step, because 1 iteration per time step results in an incorrect steady-state solution for the wake dynamics.
Compared to Wandel et al. (), Metamizer reaches a similar performance at 5 iterations per time step. At 20 or even 100 iterations, our approach is orders of magnitudes more accurate while still being reasonably fast.

#### 5.2.2 Burgers equation

The 2D coupled Burgers equation is an important non-linear PDE to study the formation of shock patterns.

$$ $F=\partial_{t}\vec{v}+\left(\vec{v}\cdot\nabla\right)\vec{v}-\mu\Delta\vec{v}=0$ (10) $$

Here, $\mu$ is a viscosity parameter. Figure 1 shows an exemplary simulation result by Metamizer that exhibits clear shock discontinuities at $\mu=0.3$. Additional qualitative results are shown in Figure 10 of Appendix A. Note that the Burgers equation was not considered during training. Mean squared residuals for different iterations per time step are shown in Table 1.

#### 5.2.3 Cloth dynamics

To describe cloth dynamics, we closely follow the approach of Santesteban et al. () and Stotko et al. (). Instead of using a classical PDE formulation, their approach relies on a spring-mass system consisting of a regular grid of vertex positions $\vec{x}$, velocities $\vec{v}$ and accelerations $\vec{a}$ that get integrated by a backward Euler scheme:

$$ $\vec{v}^{t}=\vec{v}^{t-1}+\Delta t\,\vec{a}^{t}\textrm{ and }\vec{x}^{t}=\vec{ x}^{t-1}+\Delta t\,\vec{v}^{t}$ (11) $$

To obtain the accelerations $\vec{a}^{t}$, the following physics-based loss $L$ must be minimized by Metamizer:

$$ $L=E_{\textrm{int}}(\vec{x}^{t})-\Delta t^{2}\langle\vec{F}_{\mathrm{ext}},\vec {a}^{t}\rangle+\frac{1}{2}(\Delta t)^{2}\langle\vec{a}^{t},M\vec{a}^{t}\rangle$ (12) $$

Here, $\vec{F}_{\mathrm{ext}}$ corresponds to external forces (such as gravity or wind), $M$ is the mass matrix and $E_{\textrm{int}}$ are the internal energies of the cloth specified as follows:

$$ $E_{\textrm{int}}=c_{\textrm{stiff}}\frac{1}{2}\sum_{\vec{e}\in\textrm{edges}} \left(||\vec{e}||-1\right)^{2}+c_{\textrm{shear}}\frac{1}{2}\sum_{\psi\in \textrm{right angles}}(\psi-90^{\circ})^{2}+c_{\textrm{bend}}\frac{1}{2}\sum_{ \theta\in\textrm{straight angles}}(\theta-180^{\circ})^{2}$ $$

where $c_{\textrm{stiff}},c_{\textrm{shear}},c_{\textrm{bend}}$ are parameters that describe stiffness, shearing and bending properties of the cloth.

Figure: Figure 9: Cloth with various material parameters $(c_{\textrm{stiff}}/c_{\textrm{shear}}/c_{\textrm{bend}})$ can be simulated. a) $(1000/10/0.01)$ b) $(1000/1000/10)$ c) $(1000/10/1000)$
Refer to caption: extracted/5916378/imgs/Cloth_parameters.png

Figure 9 shows results by Metamizer for various cloth parameters. The simulations were performed at 10 iterations per time step and around 185 iterations per second. Note that $L$ can take negative values and cannot be interpreted as easily as the mean squared residuals of the previous PDEs. Thus we ommited $L$ in Table 1 for our cloth simulations.

## 6 Conclusion

In this work, we introduced Metamizer, a novel neural optimizer for fast and accurate physics simulations.
To this end, we propose a scale-invariant architecture that suggests improved gradient descent update steps to speed up convergence on a physics-based loss at arbitrary scales.
Metamizer was trained without any training data but directly on the mean square residuals of PDEs such as the Laplace, advection-diffusion or Navier-Stokes equations and cloth simulations in a meta-learning manner. Remarkably, our method also generalizes to PDEs that were not covered during training such as the Poisson, wave and Burgers equation.
By choosing a proper number of iterations per timestep, a trade-off between speed and accuracy can be made. Our results demonstrate that Metamizer produces fast as well as highly accurate results across a wide range of physical systems.

In the future, further physical constrains such as Neumann or Robin boundary conditions and self-collisions for cloth-dynamics could be included. For increased flexibility regarding the underlying data structure, the network architecture could be extended to Mesh- or Graph-Neural Networks. Furthermore, multiple PDEs such as the Navier-Stokes and advection-diffusion equation could be coupled.

We believe that our approach may have applications in future numerical solvers and speed up accurate simulations in computer games or computer-aided engineering. Since gradient descent approaches are ubiquitous, our approach may also serve as an inspiration for other gradient-based applications outside of physics, such as style transfer, scene reconstruction, or optimal control.

## Appendix A Qualitative results for Burgers, advection-diffusion and Poisson Equation

In this Section we provide additional qualitative results for Burgers equation with $\mu=0.3/1$ (Figure 10), the advection-diffusion equation with $D=0.1/0.5/2/10$ (Figure 11) and the Poisson equation (Figure 12).

Figure: Figure 10: Burgers equation with $\mu=0.3/1$
Refer to caption: extracted/5916378/imgs/Burgers_results.png

Figure: Figure 11: Advection-Diffusion equation with $D=0.1/0.5/2/10$. (Domain Boundaries have Dirichlet BC = 0.)
Refer to caption: extracted/5916378/imgs/Diffusion_parameters.png

Figure: Figure 12: Solution of the Poisson equation for the electrostatic potential of two oppositely charged particles within a box with 0-boundary conditions. Note, that Metamizer was only trained on the Laplace equation ($f=0$) while, here, $f\neq 0$ at the particle positions.
Refer to caption: extracted/5916378/imgs/poisson_2_particles.png

## Appendix B Performance comparisons on a larger grid

In order to investigate the scaling behavior of Metamizer, we trained another network on a $400\times 400$ grid to solve the Laplace equation. Figure 13 shows that Metamizers advantage over traditional iterative sparse linear system solvers becomes even more significant on a $400\times 400$ grid, showing about 10 times faster convergence compared to lgmres. This is because Metamizer can be effectively parallelized on a GPU.

Figure: Figure 13: On a $400\times 400$ grid, the speed-up of Metamizer compared to other approaches becomes even more significant.
Refer to caption: extracted/5916378/imgs/laplace_400x400.png

## Appendix C Scaling parameters

In Figure 14, we show additional results how Metamizer scales $s_{i}$ for $3,5,10,20,50,100$ iterations per timestep on the Advection-Diffusion equation. For 50 and 100 iterations, a significant amount of iterations is ”wasted” to scale up $s_{i}$ for the next timestep. This could be improved by automatically scaling $s_{i}$ to a certain percentage (e.g. 20 %) of the maximum scale of the previous time step when starting with a new time step.

Figure: Figure 14: Scaling parameters for Advection-Diffusion for different numbers of iterations per time step
Refer to caption: extracted/5916378/imgs/Appendix_Scales_advection_diffusion.png

## Appendix D Laplace operator Details

Figure: Figure 15: Implementation detail for Laplace operator - example in 1D
Refer to caption: extracted/5916378/imgs/laplace_implementiation_detail.png

During our experiments, we found that naively applying the finite-difference Laplace operator results in very small gradients - even if the residuals are still quite large. Intuition on this problem is provided in Figure 15: If we want to lower the Laplacian for the center vertex (marked in black), we could increase the value of that vertex (marked in green) but we could also decrease the values of its neighboring vertices (marked in red). The same is true for all vertices, such that in the end the gradients pointing upwards get canceled out fairly exactly be the neighboring gradients pointing downwards. To avoid this issue, we detach (and thereby remove) all the neighbor-gradients (marked in red) so that we are only left with the green gradients and can make much better progress in solving the Laplace operator.

## Appendix E Video

The supplementary video demonstrates how Metamizer solves the Laplace, diffusion-advection, wave, Navier-Stokes and Burgers equations as well as cloth simulations. The time-dependent PDEs are visualized at the speed of the simulation.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 3
presentation: 2
questions: 1. The authors criticize gradient descent-like methods like AdaGrad and Adam, yet their method is based on improving gradient descent.  Can the authors speculate on how their ideas could work if based on other techniques like Adam?

2. L251-254 The authors describe how they want to use one network for all PDEs but then describe how different PDEs require different numbers of channels.  This sounds like different networks need to be trained depending on how many variables are in the problem (e.g., you'd have a 1-variable network, a 3-variable network, etc.).  Can the authors clarify this?

3. Figure 6 and the discussion around it doesn't say what equation and what boundary conditions are being solved here.  If this is Laplace/Poisson (at least without all-Neumann boundary conditions) - I don't really believe the results, because this should be a symmetric positive-definite sparse linear system that conjugate gradients should perform quite well on (and certainly better than GMRES and other slower/more general Krylov solvers).

4. With regards to the Figure 6 results, it's a little unfair to compare against constant-learning-rate variants of the other methods shown.  Of course things like SGD can't converge to 1e-33 precision with a learning rate of 0.01 - it is probably doing a very good job but the learning rate is so large that it literally cannot land on the minimum.  And it's not a mystery how to choose step sizes that enable convergence - e.g., for gradient descent, the Wolfe conditions (for instance) can inform what step size you should take at each step.

5. These results are further unfair because Metamizer - along with any other L2O methods used for these types of problems, see for instance the discussion in Kaneda et al. 2023 - are effectively acting as preconditioned solvers.  In Metamizer's case you're preconditioning gradient descent by using data about PDE losses.  So it would be fairer to compare against preconditioned CG, MINRES, GMRES, etc.  Unfortunately, the timings and iteration counts are so close that I bet if you used even just a Jacobi/diagonal preconditioner on those methods, they would all be outperforming Metamizer.

6. Figure 4 indicates that Metamizer-produced solutions may not enforce much smoothness (e.g., there is some rough noise around 20 iterations, where solutions are still not in the regime of double precision limits).  As mentioned earlier, I might venture a guess that this is due to the U-Net architecture chosen.  However, I am curious for the authors' thoughts?  This may be a useful thing to discuss in the paper, too.

Although my score leans somewhat negative (I would prefer to mark this as a 4, but it's not an option), it is mostly due to the evaluation and presentation, and I believe this is a worthwhile idea that could be published in an appropriate venue and with suitable manuscript improvements.
soundness: 2
strengths: 1. A nice advantage of Metamizer is that it's able to be applied to PDEs that weren't seen during training.  I don't know if there are other state-of-the-art methods like that, but I appreciated that feature.

2. Metamizer appears to be scale-invariant and to not need retraining for different PDEs, which again, is quite useful from a practical perspective.

3. Metamizer seems like a simple idea with a relatively straightforward architecture, which means that future research in this direction should be possible to build upon the present work.
summary: The paper proposes a meta-optimization technique for solving physics simulation problems, including linear and nonlinear PDE problems.  The proposed Metamizer framework is an iterative optimizer that uses a scale-invariant architecture to improve gradient descent updates to accelerate convergence.
weaknesses: 1. "However, to the best of our knowledge, L2O has not yet been applied in the context of physics simulations"
Unfortunately, I must inform the authors that L2O has absolutely been applied in the context of physics simulations.  Please see for example the works "A deep conjugate direction method for iteratively solving linear systems" (ICML 2023) and "A Neural-Preconditioned Poisson Solver for Mixed Dirichlet and Neumann Boundary Conditions" (ICML 2024) and the references / related work discussions in there.

2. "Data-driven Deep Learning approaches"
The authors may be interested recent works like "Toward Improving Boussinesq Flow Simulations by Learning with Compressible Flow" (PASC 2024), which show that the addition of a neural surrogate can produce results that are more accurate than what is possible with a classical numerical solver - because a solver may be using a less accurate physical model.  This is an important point to be aware of or mention, since many deep learning approaches overlook the fact that one way to improve efficiency and accuracy is to simulate things with a cheaper model and then use learning models trained on more accurate models.

2. "Unfortunately, the learned implicit representations need to be retrained for every physical simulation and thus are not real-time capable."
This criticism of PINNs does not make sense to me.  When you train a PINN, you can evaluate is very quickly.  The fact that PINNs are best when trained on a specific physics problem, does not seem correlated to the performance of evaluating (performing inference on) a neural network like a PINN.  I suggest the authors remove this criticism or clarify what they mean.

3. (Minor) In the section "Time-Integration Schemes," it is likely worth tossing in a reference to a textbook or some similar resource that discusses additional common time integration schemes like RK4 - just to avoid giving the impression that the three schemes mentioned are the only ones.

4. It is interesting the authors' technique seems to essentially just be a U-net with carefully-chosen features (Fig. 3).  With physics problems, solutions are often smooth/continuous, so practitioners often use CNNs or perhaps mix convolutional layers with U-nets.  There is no discussion of why a plain U-net architecture was chosen, and adding something like that would be helpful for readers to gain intuition into the method.

5. There's not much that can be gleaned from Figs. 7 or 8 without quantitative metrics or labeled error/scale bars.

6. L757 "This is because Metamizer can be effectively parallelized on a GPU." Sparse iterative linear solvers can be effectively parallelized on a GPU - in fact, this is one of the things GPUs are best at.  This statement makes me believe the authors are using CPU or perhaps even serial implementations of things like CG or MINRES, which is a severely (10x+) unfair comparison.

Minor:

7. L245 you call it a CNN but it's specifically a U-Net?

8. L264 You say you use gradient descent but in parentheses you say Adam - do the authors understand these are distinct numerical methods?

9. L268 "ca 6 hours" should just be "about 6 hours," or even better, just specify the time more precisely, e.g., "6 hours and 15 minutes"

10. The abstract is a little long and split into a few paragraphs; I would recommend condensing things a bit.

11. Formatting mistakes like not doing LaTeX quotes properly - the authors should know better than this in a venue like ICLR

12. L235 citation not formatted right (\citet not \citep)

13. L277-281 (for instance - there are other places for this) Put commas after e.g., this is a grammar thing

14. L434 probably a typo in the table here

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: ## Missing details/questions

**Q1**: Insufficient details on residual computation (148-154, 175), specifically handling spatial resolution and Navier-Stokes pressure. More information should be provided in the appendix (related to **MW3**).

**Q2**: In (198-211) you say that the main problem with directly optimizing the gradient is the scaling, therefore you introduce a scale-invariant optimizer. What about local minima? The solution of the physical problem is the global minimum of the problem, so getting stuck in a local minimum is a potential issue that is not necessarily solved by scale invariance.

**Q3**: Is there no cold start issue? In the **Training** paragraph (256-269) it is stated that the training pool is initialized with randomized solutions. How exactly?. Fig. 6 shows that the gradient-based optimizers (Adam, SGD, etc) fail to converge. If the gradient signal is so poor how can the UNet learn from it? How does this relate to the scale vs local minima issues from Q1? Can Metamizer, once trained, minimize for randomized solutions drawn from a different distribution than the one it was trained on?

**Q4**: I appreciate the wall time and GPU usage training information (269). How many update steps were performed on the model? How many data samples did it see (presumably the number of update steps * batch size)?

**Q5** I find the description of the time integration and how it is used in the context of Metamizer unclear. Is the integration scheme used to define the discretization used of the time derivative? E.g. for forward Euler it would be approximated as $u_t' = \frac{u_{t+1} - u_{t}}{\Delta t}$ 

**Q6** Why is time integration needed at all? Can the residual computed on the full trajectory not be minimized for all the states at once?


## Things to improve the paper that did **not** impact the score:

Abstracts usually have a single paragraph.

Space before comma on line 93.


## Refrences

1) Brenner, S. C., & Scott, L. R. (2008). The Mathematical Theory of Finite Element Methods. In Texts in Applied Mathematics. Springer New York. https://doi.org/10.1007/978-0-387-75934-0 
2) Chassaing, P. (2022). Fundamentals of Fluid Mechanics. Springer International Publishing. https://doi.org/10.1007/978-3-031-10086-4


## Edit history
1. Confidence reduced from 4 to 3 in light of the discussion with the authors.
soundness: 1
strengths: The problem the paper tries to address is of interest and has wide applicability.
summary: The paper presents a neural-network-based meta-optimizer, Metamizer, designed to minimize the residuals of physical systems, especially partial differential equations (PDEs), in pixel space. The authors conduct experiments comparing the performance of Metamizer to both gradient-base optimzers and linear system solvers. Metamizer is evaluated on different PDEs and the results presented.
weaknesses: The paper has 3 main weaknesses:
- **MW1**: Lack of quantitative evaluation. Only qualitative images and residuals are presented, making it difficult to assess the method's performance objectively. This issue affects all the experiments.
- **MW2**:  No baseline comparisons with other Learning to Optimize (L2O) methods, classical numerical solvers, or neural PDE methods.
- **MW3**:  Limited experimental scope. Experiments use only fixed-size grids, 100x100 (331) and 400x400 (775), which restricts generalizability and does not test varying resolutions, where new dynamics may arise, especially in the Navier-Stokes equation.


## Main weaknesses

The proposed method can be viewed as a meta-learning method, a PDE solver, or an emulator but it lacks comparisons with other methods from these fields. Ideally, the method would be compared against methods from each of these fields. I acknowledge that such an exhaustive comparison involves a considerable effort, but the paper omits even a basic evaluation against standard methods in these areas.

The authors assert that L2O has not been applied to physics (99) but fail to demonstrate how Metamizer performs in the physics context against other L2O methods. Furthermore, the paper’s accuracy claims lack quantitative evidence, as it provides only qualitative images and residuals. This is inadequate for judging solution quality, especially when standard neural PDE literature benchmarks against ground truth obtained from numerical methods. Moreover, the lack of comparison to classical numerical methods is striking, since a comparison is performed against linear solvers for a linear system built in pixel space. Effectively, the paper is comparing against different linear solvers for the finite element method.

In short, the experimental evidence does not convincingly support the method's efficacy or comparative performance and extensive existing literature has not been engaged with.

## Additional Weaknesses

**W1**: The claim that no hyperparameters are used (68-69) is inaccurate, as $s_0$ is a hyperparameter (238).

**W2**: The claim on (332-333) is misleading. In the deep learning literature, accuracy typically reflects error metrics like MSE between the ground truth and the prediction. In this paper, accuracy has a different meaning (value of the residual) so it is not directly comparable. This is related to **MW1-2**.

**W2**: The accuracy claim on (332-333) is misleading; in deep learning, accuracy typically reflects error metrics like MSE. It is not possible to directly compare it to the claimed MSR.

**W3**: Misleading claim about linear solvers and their usage. In Fig. 6 Metamizer is compared to linear solvers, and in (327) it is claimed that linear solvers are limited to linear PDEs. This is only true in pixel space. Most PDE-solving techniques, such as the FEM, rely on solving a linear system derived from the equation through linearization [1]. Related to **MW2**.

**W4**: The caption of Fig. 8 claims that Metamizer can simulate a wide range of Reynolds numbers and shows snapshots for $Re=1$, $Re=50$ and $Re=200$. This range is insufficient compared to the typical Reynolds number range (1 to millions) in fluid dynamics, missing the complex behaviors at higher Reynolds numbers [2].

**W5**: No limitations whatsoever are discussed or acknowledged.


-------
Given all of the above, I recommend the paper be rejected.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: Closely related to the topics above, can the authors shed light on the following topics?

- How does the approach differ from the scale-invariant learning proposed in the NeurIPS'22 paper above?

- Does the theory there affect or relate to the current submission? 

- How are the scaling factors s_i computed, and how accurately are they inferred?

- Have the authors considered Newton methods, that approximate the inverse Hessian for "scale invariance"? I think this would also be worth a discussion, at least.

- The paper focuses on solving optimization problems in terms of PDE states. How would this extend to a second neural network, i.e. backpropagating further than just the u?
soundness: 3
strengths: The central goal of addressing "scale invariance" is not new, but the proposed learning of scale invariance is new if I understood the approach correctly.

The breadth of results is also great to see.
summary: This paper targets an optimizer for training deep neural networks for improving the solving of tasks related to PDE solving. The goal if scale invariance is identified as a central challenge for physics-related learning tasks. The paper then proposes a U-net that is employed to infer the scaling of an update in order to improve conditioning of the neural network update.

The importance of scale invariance has been studied before for physics problems, e.g. in Holl et al., Scale-invariant learning by physics inversion, NeurIPS 2022. I just double checked , and e..g fig 1 there and figure 2 here are actually have a lot in common. This previous paper in the end takes a different approach, so I think it's no direction competition, but the unclear parts of the current submission definitely lead to a few questions.

The paper evaluates this idea for several canonical PDEs (Poisson & Laplace, advection-diffusion), as well as more complicated ones (waves, Navier-Stokes). Interestingly, cloth dynamics are included as well. While the rang of experiments is very nice, the motivation in terms of theory behind the approach remains unclear. Nonetheless, I think it’s an interesting and novel direction. So I liked the paper and its direction, but several unclear topics remain, which I hope the authors can address during the rebuttal phase.
weaknesses: One central part that was unclear to me was how the scale factors for learning are actually precomputed. Is this simply a scaling of the reference solutions from the solver? As this is the central component that is supposed to make the method work, I think it would be important to be clear where this comes from. Correspondingly, I was missing an evaluation in terms of how well the network actually learned to predict these scaling factors. The results seem to only show the loss in terms of the state.

The Holl'22 paper above also provides some theory on how updates from scale invariant solvers improve learning. Unfortunately the submission here does not take up on these topics, and focuses on experiments.

The appendix is also very brief on how the different learning tasks were actually set up.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 3
questions: - Does CG fail to converge in Fig. 6? Is the matrix not symmetric positive definite?
soundness: 2
strengths: - The idea of applying meta-optimizers to learning-based PDE solvers looks new to me. I think this is an interesting thought worth further discussion and exploration.
- The technical method is simple, straightforward, and easy to implement.
- The paper writing is good; it articulates its technical details very well.
- The paper considers a diverse set of PDE problems, including N-S equations and 3D cloth simulation (albeit with limited features).
summary: The paper studies solving partial differential equations (PDEs) numerically. The proposed framework builds upon a finite-difference spatial discretization of the unknown fields on a regular grid. Next, it minimizes a residual loss similar to that in the physics-informed neural network (PINN) via gradient-based optimization. The new idea comes in a neural network that predicts the update rule (scaling factor and direction) used in each step of the gradient-based optimization. The paper demonstrates this idea on several linear and nonlinear PDE problems.
weaknesses: **Introduction/Related Work**
- The paper claims numerical methods for solving PDEs “are computationally
expensive”/“are still relatively expensive to compute”/“require high computational resources.” This argument is generally valid for large-scale problems. However, for the small-sized problems demonstrated in this work (400 x 400), Chances are that numerical methods are not as expensive as the paper indicates and may be faster than Metamizer (See Poisson/Laplacian problem below). I suggest the paper evaluate its method on problems of much larger sizes, e.g., on >=1024^2 (2D) or >=64^3 (3D) grids.

**Foundations**
- Equations 3 and 4 are not new and should cite the original PINN paper. It is probably also worth mentioning that minimizing Equation 3 may get stuck in a local minimum (grad L = 0, but L is nonzero), which is not a solution to the original PDE. For several problems (linear second-order PDEs) demonstrated in the experiments, a variational loss from the PDE’s weak form is a standard and a naturally better candidate for the loss design. The incremental potential in Eqn. 12 is one such example that highlights this point.
- The issue in Fig. 2 and the main intuition for taking gradients from previous steps into account have been extensively studied in numerical optimization methods. The former is typically handled with (exact) line searches, and the latter is discussed in (nonlinear) CG, L-BFGS, etc. (“Numerical Optimization” by Nocedal & Wright.) It is still very interesting to contemplate the possibility of training a neural network to exploit gradients from previous steps, but the paper should discuss its relation to these classic numerical optimization methods and compare them in an experiment.

**Experiments**
- The Poisson/Laplacian example has the most complete analysis, but the main results in Fig. 6 look dubious for several reasons:
1) It looks like Metamizer has used an Nvidia RTX 4090 GPU (line 330 and Appendix B), so it should be compared with linear solvers implemented on the same GPU. However, I don’t recall scipy having a GPU implementation of them. Please correct me if I am wrong.
2) The state-of-the-art GPU implementation of a Krylov subspace solver (CG, GMRES, etc.) is likely to be better than both the scipy baselines and Metamizer reported in Fig. 6 and Appendix B. Back in 2005, deploying a standard CG algorithm on a GPU could solve wave/Poisson equations on a 1024^2 grid at a rate of 10-30 fps (Chapter 44 in “GPU Gems 2” edited by Pharr & Fernando and the referenced papers), already arguably faster than the speed of Metamizer reported in the much smaller 400^2 problem. I feel this paper didn’t identify the right baseline for comparison.

- The wave equation example does not compare Metamizer’s performance with a baseline. I suggest comparing it with the CG-GPU baseline mentioned above.

- Eqn. 12 in the cloth example is not new and should cite its original source, e.g., Stotko 2024. This equation can be traced back to at least “Optimization Integrator for Large Time Steps” by Gast et al. in 2015, if not earlier.

- “Note that L can take negative values and cannot be interpreted as easily as the mean squared residuals of the previous PDEs. Thus we omitted L in Table 1 for our cloth simulations.” In this case, grad L = 0 is the right indicator for successfully solving the PDE. There should be a table reporting grad L after certain iterations.

- On the algorithmic side, I feel some classic approaches for exploiting gradients from past iterations should be considered for comparison to establish the advantage of the network module in Fig. 3a. Two candidates are nonlinear CG and L-BFGS (with history size = 2). These methods do not require training and can be directly used to minimize Eqns. 3 and 4.

Meta-review: **Summary and strengths** The paper presents a learned solver for partial-differential equations (PDEs) which minimizes a physics-based loss function. The method achieves high accuracy on a wide range of PDEs, such as Laplace, advection-diffusion and incompressible Navier-Stokes equation. Unlike other existing neural solvers, at test time, Metamizer can be applied to solve new PDEs not seen during training. such as Poisson, wave and Burgers equation. The approach does not require having a dataset of PDE solutions, as the approach uses the solutions from the previous iterations as the training data.

**Weaknesses** Reviewers pointed out insufficient comparison to the classic solvers like CG, other GPU-based solvers, comparison to other Learning-to-Optimize (L2O) approaches for PDEs. During the rebuttal, the authors added extensive comparisons to address these points.

**Decision** The decision is Acceptance (poster). The paper proposes a simple, yet universal technique that can solve a wide range of PDEs with a single pretrained model, including unseen PDEs, having competitive performance to classic solvers. Reviewers shqZ, EsdF, 8DMR agree that the approach is valuable for the ICLR community and can handle various PDE problems.

**Final version** We suggest the authors to (1) add a specific comparison between a GPU Metamizer to GPU PCG comparison in fig. 17. (2) make more precise claims about the advantages of the method

Decision: Accept (Poster)