Riemann Tensor Neural Networks: Learning Conservative Systems with Physics-Constrained Networks
Authors: Anas Jnini
Lorenzo Breschi
Flavio Vella
Venue: ICML 2025 poster (2025)
Keywords: Divergence-free symmetric tensors
Riemann Tensor Neural Networks
neural PDE surrogates
physics-constrained neural networks
PDEs
inductive bias
fluid mechanics
Primary area: applications->chemistry_physics_and_earth_sciences
TL;DR: We propose RTNNs, neural architectures that parameterize divergence-free symmetric tensors (DFSTs), encoding conservation laws like mass and momentum at the architectural level to model conservative dynamical systems .

Abstract: Divergence-free symmetric tensors (DFSTs) are fundamental in continuum mechanics, encoding conservation laws such as mass and momentum conservation. We introduce Riemann Tensor Neural Networks (RTNNs), a novel neural architecture that inherently satisfies the DFST condition to machine precision, providing a strong inductive bias for enforcing these conservation laws. We prove that RTNNs can approximate any sufficiently smooth DFST with arbitrary precision and demonstrate their effectiveness as surrogates for conservative PDEs, achieving improved accuracy across benchmarks. This work is the first to use DFSTs as an inductive bias in neural PDE surrogates and to explicitly enforce the conservation of both mass and momentum within a physics-constrained neural architecture.

Riemann Tensor Neural Networks: Learning Conservative Systems with Plysiae-Constrained Networks
Anas Jnini"
Lorenzo Breschi Flavio Vella
Abstract
Physical Inductive Biases in Machine Learning. A cen- tral limitation of generic neural models is their lack of built- in physical intuition. While convolutional or attention based layers successfully exploit certain data symmetries (e.g.. translation invariance), they do not automatically enforce fundamental physics, such mass conservation energy preservation. Physics-Iffrmmed Informed Neural Networks (PINNs)
Divergence-free free symmetric tensors (DFSTs) are fundamental in continuum mechanics, encoding conservation laws such as mass and momentum conservation. We introduce Riemann Tensor Neu- ral Networks (RTNNs), novel neural architec- ture that inherently satisfies the DFST condition to machine precision, providing a strong induc- tive bias for enforcing these conservation laws. We prove that RTNNs can approximate any suffi- ciently smooth DFST with arbitrary precision and demonstrate their effectiveness surrogates for conservative PDEs, achieving improved accuracy across benchmarks. This work is the first to use DFSTs as an inductive bias in neural PDE surro- gates and to explicitly enforce the conservation of both mass and momentum within physics- constrained neural architecture.
(Lagaris et al., 1998; Raissi et al., 2019h; Cai et al., 2021; Haghighat et al., 2021; Hu et al., 2023) address this gap by adding PDE residuals and boundary conditions as soft con- straints in the loss function. PINNs have been successfully deployed many PDE problems (Karniadakis et al., 2021; Jnini et al., 2024b), but the "soft penalty" approach can lead to suboptimal enforcement of conservation laws and stiff optimization (Wang et al., 2021). Consequently, there is growing interest in hard or explicit constraints that guar- antee PDE structure a priori (Richter-Powell et al., 2022; Greydanus et al., 2019b; Cranmer al., 2020; Jnini et al., 2024a; Liu et al., 2024).
Divergence-Free Symmetric Tensors in Physics and Mathematics. Divergence-free symmetric tensors (DF- STs) are a special class of tensor fields characterized by vanishing row-wise divergence and inherent symmetry in their indices. In (n + dimdennsio space-time do- main, where represents the spatial dimensions and the additional dimension accounts for time, these tensors fre- quently arise as stress or momentum flux tensors in various fields, including fluid dynamics, elasticity, kinetic theory, and relativistic hydrodynamics (Serre, 2018). For instance, the Navier-Stokes stress tensor or the compressible flux matrix can be represented an (n 1) x (n + 1) tensor S satisfying V S = 0. This formulation unifies the conserva- tion of mass and momentum under single flux-divergence constraint. Despite their natural alignment with PDE-based conservation laws, existing studies have not explored lever- aging DFST structures directly as a neural inductive bias within neural network architectures.
1. Introduction
Partial Differential Equations (PDEs) Partial Differen- tial Equations (PDEs) are central to the mathematical model- ing of complex physical systems, including fluid dynamics, thermodynamics, and material sciences. Traditional nu- merical methods, such as finite element and spectral meth- ods, often require fine discretization of the physical domain to achieve high accuracy. These approaches become computationally expensive, particularly in engineering ap- plications where systems must be solved repeatedly under varying parameters or initial conditions. Recent advances in machine learning (ML.) have shown promise in addressing these challenges by leveraging neural networks (NNs) as po- tential alternatives enhancements traditional numerical solvers (Kovachki et al., 2021; Li et al., 2020).
Department of Information Engineering and Computer Sci- ence, University of Trento, Trento, Italy. Correspondence to: Anas Jnini <anas <anas,jninie jnini @ unitn.it>..
Our Contributions. In this work, we present a novel approach to embedding fundamental conservation laws di- rectly into neural network architectures through DFSTs Our primary contributions are follows:
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constraineed Networks
1. Architectural Design of Riemann Tensor Neural Networks (RTNNs): We introduce RTNNS, class of neural architectures specifically designed to gener- ate DFST fields. RTNNs are tailored for approximating individual DFSTs, ensuring the divergence- condi- tion, V S = 0, is satisfied to machine precision.
ensures strict conservation and symmetry. To the best of knowledge, this work is the first to use enforce the conserva- tion of both mass and momentum at the architectural level for surrogate modeling.
2. Background and Theory
2. Theoretical Guarantees for RTNNs:
We establish theoretical foundations for RTNNs by proving their uni- versal approximation capabilities for any sufficiently smooth DFST.
Notation (Preliminaries) Let n denote the spatial dimen- sion and ? IR the spatial domain. The space-time do- main is f2 [0. x Ω C R=*1, where E [0, T) and x E Ω. For a function [(I, x), the augmented gradient is f f = (0rf. Pe..,...,0r Gx,,f) f), and the spatial gradient is
3. Empirical Validation and Comparative Analysis:
Vf fo........00
We reformulate several benchmark problems within the DFST framework and conduct numerical experi- ments that demonstrate RTNNs consistently improve performance of PINNs in accuracy when used as sur- rogate models for conservative PDEs.
Divergence- Freee Symmetric Tensors in Continuum Me- chanics. We begin by introducing the class of divergence- free symmetric tensors, which encode either the conser vation of mass and momentum in classical mechanics or energy and momentum in special relativity. A tensor field
In the following sections, we review DFST-based PDE for- mulations, describe the proposed neural architectures, and present experimental validations on benchmark problems.
S: Ω2 - [(++1)x(++))
is said to be symmetric if S and divergence-free if it satisfies
1.1. Related works
Divergence-Free Symmetric Tensors in Mathematical
Div (S): desio + > Ong = 0, E nt..
Physics. A large body of work by (Serre, 1997; 2018; 2019; 2021) has established the fundamental importance of DFSTs in continuum mechanics and kinetic theory. These tensors encode conservation principles for mass and mo- mentum (in classical fluid dynamics) or (in relativistic hydrodynamics), and are present in models ranging from Euler or Boltzmann equations to mean-field (Vlasov-Poisson) descriptions of plasmas and galaxies. Al- though DFSTs have been investigated in PDE theory, prior investigations have largely focused on analytic or qualitative properties. To the best of our knowledge, no existing work leverages DFSTs explicitly as a numerical method or as an architectural inductive bias in machine learning frameworks.
Many classical PDE systems, including compressible or in- compressible fluid flow, elasticity, and shallow- water mod- els, can be expressed in this flux-divergence form by suitably choosing S. canonical representation in fluid mechanies is:
m m & m + σ
S
(1)
P
where Or + R denotes the mass density, m 22 + Rn the linear momentum field, and a : 2r -> T2 the stress tensor. Enforcing (S) = 0 then unifies mass and momentum conservation, while additional constraints (e.g., constitutive laws or energy equations) can further spec- ify couple p and m.
Hard-Constraints in Scientific Machine Learning. Be- yond the classical physics-informed approach of adding PDE residuals soft constraints in the loss (Raissi et al., 2019a; Karniadakis et al., 2021), there is growing inter- est in incorporating hard constraints or specialized struc- tures into neural networks. For instance, (Hendriks et al., 2020) investigate linearly constrained networks, (Richter- Powell et al., 2022) impose iantinuity-equaiion constraints via divergence-Hig
Motivation for Neural Parametrization. Although one can penalize the residuals of the condition Div = in a soft-constraint manner (e.g., through terms in the loss func- tion), this approach does not guarantee the satisfaction of the divergence-frr condition, especially when optimization is challenging regularization terms underweighted; this is often the case in the non-linear regime that we are considering, as described in (Bonfanti et al., 2024).
vector fields, and several recent meth- ods aim preserve energy momentum (Greydanus et al., 2019a; Cranmer et al., 2020). These efforts reflect a broader push in machine learning to embed domain-specific priors, thereby improving stability and generalization (LeCun et al., 1998; Giles& & Maxwell, 1987). Our work similarly encodes the PDE structure "at the network level" via DFST, which
Instead, we propose embedding the divergence-frr ffre prop- erty directly into the neural network architecture. A pri- mary motivation for this architectural choice, beyond the
2

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Plysies-Constrainee Networks
immediate benefit of ensuring strict conservation, more fundamentally address how coupled physical quantities are represented. In many physical systems, such as those in fluid dynamics, variables like density, velocity, and pressure are not independent but are deeply interrelated. Conventional neural network models that treat these as separate output channels can struggle to capture these intrinsic physical and geometric correlations for instance, the interdependence between velocity components the unified response of density and velocity to pressure gradients. By designing our network to output DFST, we aim to inherently model these quantities as components of a single, unified field, linked directly by the underlying conservation laws. Consequently, such hard-coded constraint ensures strict conservation (c.g.. of mass and momentum) to machine precision, while also providing better physical consistency and stronger inductive bias.
The space of divergence-free symmetric tensors Sab on the flat manifold is the image of the map:
Sab V. K Kachd
(6)
where Kabed is any (0, 4)-tensor satisfying the symmetries (2)-(4) Any Sab can be expressed as:
Sab Σ
V Ci::
ISiSism
where cij are smooth scalar functions.
Proof. The proof is presented in Appendix A.
2.2. Riemann Tensor Neural Network
In the setting where we wish to approximate a divergence- free symmetric tensor field
Proposed Approach. The following sections introduce nettlllntworr-bbseeee construction that guarantee the output is frrergence-free symmetric tensor. We prove that our approach can approximate any sufficiently smooth DFST to arbitrary accuracy, thus offering a robust way to integrate conservation principles into neural PDE solvers.
S: - Raxx
(8)
we define a Riemann Tensor Neural Network as follows.
Definition 2.2 (Riemann Tensor Neural Network'). Sup- pose:
2.1. Constructing Divergence-Free Symmetric Tensors on a Flat Manifold (DFSTs)
Ω C R" is flat domain (manifold without boundary or with suitable BCs).
Theorem 2.1 (Representation of Divergence-Fres Symmet- ric Tensors on a Flat Manifold). Let V be an n-dimensional real vector space with a fixed basis {ea Ja-i. and let {e-)== denote the corresponding dual basis of V*. Let A²V* denote the space of 2-forms on V. Consider the space of all (0, 4)-tensors Koted defined on a flat manifold equipped with a
(T((t)) } is a finite non-traine non-trainable basis of Riemann- (0, 4)-tensors as defined in Theorem 2.1,
m(m+1)
NNo IR* + R 2 is a multilayer perceptron (MLP) with tifeeeliffeeeetiabble activations, whose in- put r E IR* indexes points in Ω.
connection V, satisfying the
following symmetries:
Then an RTNN for the single-field case is constructed by:
1. Antisymmetry within index pairs:
1. Scalar coefficients {c: 0)): The MLP NNg out- puts m(m+1) sealar functions { (Cay(z )sssssmm
K(oblydd 0,
(2)
2
K Kab(ad) = O.
(3)
2. Hessian Computation: For each cij (x; (z; 0), compute the Hessian components 8 add (z)(2; 0) via automatic differentiation.
2. Symmetry between pairs:
3. Tensor Field Construction: Define the (0, 2 -tessor field
Kated = Kedah-
(4)
Let (wi, (wi,...,Wim, Wm be a fixed basis of A2V", where m m(n-1) Define the tensors:
So(x)
Σ 1siSi Sm
((t)) achd ar (x)(2; 0).
(9)
2
((t.d) == wilc". A etwyy A A + wey A wille A ea), (5)
By design, So(x) is row- divergence- and sym- metric for all x S.
where e Aeg is the wedge product of the dual basis elements e,, and eg.
We call Sg Riemann Tensor Neural Network. It provides parametric approximation So to single DFST on Ω
3

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Remark 2.3. Although the network outputs scale as m(m--
ery involves dividing by p. leading to instability when p is initialized around small value (Richter-Powell al., 2022). To address this, we add an identity matrix to Sp, ensuring is initialized near 1 without violating DFST con- straints. Throughout the experiments in this paper, we use the Least-memory BFGS (Nocedal & Wright, 1999) opti- mizer due to the highly non-linear nature of our problems. Additionally. LBFGS efficiently approximates second-ordes curvature information, facilitating effective optimization in the complex, non-linear loss landscapes encountered in our experiments. The challenges of optimizing in such regimes and the importance of second-order optimizers have been well documented in the literature (Jnini et al., 2024b; Müller & Zeinhofer, 2024; Bonfanti et al., 2024).
1), the problem setup is generally overparameterized. De- pending the application, certain scalar functions can be set to zero without violating the DSFT condition, provided that the number of basis functions exceeds the degrees of freedom.
2.3. Universal Approximation Theorem for RTNN
Theorem 2.4 (Universal Approximation for RTNN), Lel Q C. IR" be a bounded (hence compact) domain, and sup- pose S: Ω + [Rax" is a C2 -mmoohh divergence-free, sym- metric tensor field on Ω2 Then for any E >> 0, there exists a Riemann Tensor Neural Network Sg such that
sup FE
[[5(z)) Sa(z)]]Fro < <,
(10)
Code implementation and public repository Our code has been implemented using the JAX library (Bradbury et al., 2018). Our implementation is publicly available at https: /
where ||. ||Fro denotes the Frobenius norm on Bux, In particular. Sy also remains divergence-free and symmetric on Ω
Riemann- Teeinn--enoo---eee Neura Networks.
Proof. We present the proof in Appendix A.3
3.2. Pedagogic Example: 2D Isentropic Euler Vortex
3. Methodology and Applications
For this pedagogic example, we simulate a 2D isentropic Euler vortex a smooth, rotational flow solution to the Euler equations that accurately captures vortex dynamics and is commonly used as a standard benchmark for evaluating the accuracy of numerical solvers-over a spatial domain Ω [0, Lg] x [0, Ly] and time interval |o, T). We employ well-defined analytical initial and boundary conditions that are detailed in Appendix B.1.
The preceding sections established the theoretical foun- dation of divergence-free symmetric tensors DFSTs and RTNNs as a rigorous approach enforeing conservation laws within neural architectures. In the following sec- tion, we move from theory to application, showcasing how RTNNs can be employed to model various conservative sys- tems, including the Euler and Navier-Stoke, equations, and Slggneto-YydrodyyamiisMMM
Let > 0 denote the density, (u, v) the velocity field, and the pressure. For isentropic flow with y > 1, the pressure is given by p(p) = kp". The governing equations, including the energy equation, are:
3.1. Efficient Implementation and Practical Considerations
We employ Taylor-Mode Au- tomatic Differentiation, which propagates Taylor coeffi- cients through the network by treating the computational graph as an augmented network with weight sharing. This approach effectively reduces redundant computations as- sociated with higher-order derivatives, significantly accel- crating the training process of RTNNs. Additionally, for in Section 3.4, we utilize Separa- ble Physics- Physics-Informed Neural Networks (SPINNs). SPINNs decompose PDE residuals into per-axis evaluations, facili- tating efficient differential operations on large-scale regular grids (Cho et al., 2023).
Automatic Differentiation
Governing Equations The 2D compressible Euler equa- tions are:
0.p + V (pu, = 0,
(11)
(.p))) + V (pu? .puv) -
(12)
0,(p)) + V. (puv, poo) = -OyP,
(13)
0, + V ((E -p)u, (F + p)u) = 0,
(14)
where E 7-1 + p22 2 + u2) is the total energy.
DFST Formulation and Decomposition of a. Rewriting (11)~(13 in flux-divergence form, we express the system as:
Optimization and Stability Our method models densi- ties and momenta instead of velocity fields, velocity recov-
V GZz S =0 0, where
pu [ou)² + Tza (pu)
pv (ou)pp))
So named because the underlying (0, 1) 1)-tensors share index symmetries with the Riemann curvature tensor in di differential ge- ometry.
S = pu
(pa) P
pe
+ dyy
4

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Here, the stress tensor is:
While the parametrized fields exactly satisfy the DSFT con- straints, they are only solutions the momentum equations if the stress tensor satisfies the zero deviatoric constraints, which we can penalize in the loss function in addition to the boundary and initial terms.
deviatoric part
where = p(p) represents the isotropic pressure contri- bution. The absence of a deviatoric term reflects the as- sumption of inviscid flow. The divergence-free condition
We clarify the equivalence of the zero deviatoric constraint with the momentum residual in Appendix Section A.4.
V = 0 enforces:
Loss Function. To train the RTNN and ensure that the modeled tensor So adheres to the governing equations and boundary conditions, we define an objective function:
mass conservation, and
momentum conservation in and y-directions.
(()) LBC + Lic + Lo + Le.
Additionally, the energy equation is given separately as:
where LBC penalizes deviations from the prescribed bound- ary conditions while Lıe enforces consistency with initial conditions. The term C. ensures that the stress tensor re- mains purely isotropic by penalizing the magnitude of the deviatoric component, (og"x112. it's inclusion is equivalent to penalizing the momentum equation residual as shown in in A.4. Finally, LE minimizes the residual of the energy equation, measured ||0,Eg + V ((Eg + Po) up)) l? ensuring that the total energy properly conserved within the system. All loss terms are formulated in the least squares sense.
0 . (E + +))) 0,
where the total energy E' is:
== 7-1
+ J0lu² tz))
RTNN Parametrization. To approximate solutions of (11)-(14), we define family of tensors So using RTNNs The parametrization proceeds follows:
1. RTNN Parametrization of Sg: Let Sg denote an RTNN as described in Section 2.2. By construction, Sy is divergence-free and symmetric in Qr.
Experimental Setup. For the neural network training, we sample 500 interior collocation points within Ω x [0, T) to enforce the residual constraints of the governing PDE. Ad- ditionally, 100 boundary and initial condition points are sampled to impose the prescribed constraints. Our RTNN model is parameterized by Multilayer Perceptron (MLP) with 4 hidden layers, each containing 50 neurons. Train- ing is performed entirely without labeled data, model is validated against the analytical solution of the isentropic Euler vortex evaluate accuracy.
2.
Extracting Physical Fields: We can interpret Sg in block form. From it, we read off::
Pe = (So)o.0.
(ppl)) = (Sa 1:2.0.
σg = 1S,,1 (pous, )ove) (paus, pove)
3. Zero-Deviatoric
Constraint and Energy Parametrization: We can parametrize the pres- sure by decomposing the stress tensor into isotropic and deviatoric parts:
We benchmark RTNN against two methods: (1) the standard PINN approach and (2) Neural Conservation Laws (NCL) that enforces exact mass conservation (Richter-Powaareee et al., 2022). Both methods use similar MLP architectures for fairness. Performance is evaluated terms of median av- erage relative L error on all fields and simulation time, We train all three methods using 200,000 iterations of the L-BFGS optimizer. We follow the loss scheme presented in this section, while training both PINN and NCL using PDE residuals penalized in the loss.
OB IoI + 0
where:
1 pe
2 to)))
odev oe Pol.
(15)
Additionally, we parametrize the energy as:
Results and Discussion. Table 1 summarizes the results, while Figure presents the evolution of the relative 7,2 error over simulation time. Our RTNN significantly outperforms both PINN and NCL, achieving a median relative I er- ror that two orders of magnitude lower than PINN and four orders of magnitude lower than NCL. Furthermore,
Eg
Pe
po - + vz)).
7
For an inviscid isentoopcc vortex, we enforce the con- straint ode = 0 during training.
5

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physics-Constrainct Networks
Method
Relative L2 Error
Wall Time (s)
while, equation (17) balances the convective, pressure, and viscous forces within the fluid. The term vÂu specifically models the internal fluid friction due to viscosity.
RTNN
9.92e-05
596.34
NCL
3.87c-01
2008.57
PINN
3.82e-02
365.67
DEST Formulation and Stress Decomposition. These equations can be expressed in a divergence- symmetric tensor (DFST) form. We define a tensor
Table Comparison methods the Euler experiment, report- ing the median relative L2 error and median wall time across five independent training runs with different random seeds.
10
S =
uu
10-4
such that
10:7
V = 0.
10:1
RUNN
Here, the term u @ u represents the convective flux, while the total stress tensor decomposed as
MINN
10
50
100
130
200
250
300
350
pl dav
Simulation Time (5)
RTNN (L=0.5)
PINN (t=l.5)
0.016
0,045
where the deviatoric part o devv captures viscous stresses via
0.014
0.040
Enor Magnitudi
0.012
0,035
0.010
0.025
0.030
dev = vvvu++(VVi (Vu)))
0.020
0.006
0.015
0.004
0.010
0.000
0.000
0,005
Exact incompressibilityy To enforce exact incompress- ibility (i.e., Soo = 1). observe that any contributions to So0 come specifically from basis 2-forms containing ej- Con- sequently, by choosing the corresponding coefficients oj to vanish whenever the wedge product includes eğ in Tabeds we ensure that Soo consists only of the identity term we added for stability, thus achieving exact incompressibility.
Figure Comparison training dynamics and error fields RTNN, NCL, and PINN for the Euler experiment. Top: relative T.2 error evolution over simulation time, showing median (lines) and IQR (shaded). Bottom: Error fields for RTNN (left) and PINN (right) 0.5
Following Section 3.2, we em- play an RTNN Sş to represent the solution. From its block structure, we extract the physical fields as:
RTNN Parametrization.
RTNN demonstrates stable convergence while maintaining competitive training times.
Pe 1, (us, Uo, wp,we) = (Sg) 1:n,0>
3.3. Incompressible Navier-Stokes Equation
We next consider the incompressible Navier-Stokes equa- tions, which include additional viscous forces compared to inviscid Euler flows.
ao (So) (uo: up, wo) & (uo, Vg, wa).
Here, Po represents the density field, while (uo, Vos wo) cor- respond to the velocity components. The viscous stress tensor D3 can be computed via automatic differentiation applied directly to the velocity field. Instead, we define:
Governing Equations. The incompressible Navier- Stokes system in the n-dimensional case given by:
Doey = EVuu + (Vus)T).
V u-= 0, 0,u + (u V)u = - vAu,
(16)
(17)
Although this introduces additional computational complex- ity, it can be mitigated using Taylor-mode automatic dif- ferentiation. The stress tensor can be decomposed into an isotropic part and a deviatoric part:
where u = (u, U, w) represents the velocity field, denotes the pressure field, and > 0 is the kinematic viscosity. Equation (16) enforces the incompressibility condition, en- suring that the divergence of the velocity field zero. Mean-
= PoI Do
-tr(ae).
6

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Reference Solution
Error (RTNN)
Error (NCL)
2
1.0
0
-2
2
2
4
Figure 2. Qualitative comparison of predicted velocity fields for the unsteady Cylinder Wake flow Rt 100 and representative time instance = 7.0. Left: Reference solution showing the reference velocity magnitude ||]u||]. Middle: Absolute error the RTNN prediction, e |Uueer URTNN||2. 2/2. Right: Absolute of the NCL prediction, ||uRer usct. 2. This visualization clearly demonstrates RTNN's superior accuracy in capturing the complex vortex shedding structures with significantly lower emor compared to NCL.
Viscous Residual and Loss Function. To enforce mo- mentum balance, we focus on matching the RTNN-derived deviatoric stress ogev with the velocity-based viscous stress v( (uu + (Vua) Define:
around NACA 0012 airfoil. The steady-state problem is addressed by treating time as a dummy dimension set to zero in the forward pass. The computational domain is discretized with 40,000 collocation points enforce PDE residuals and boundary conditions, alongside 2,000 labeled data points obtained from an in-house Averagee Navier-Stokes (RANS) solver to supervise the training. An MLP consisting of 4 hidden layers and 50 neurons per layer with Tanh activation functions is employed. The model undergoes 50,000 iterations of the L-BFGS optimizer. Validation conducted on the 14,000 collocation points eval- uate accuracy. Further setup details and visual results provided in Appendix B.4.
Rigivisoous = dev - DDo
Because pressure can act as a scalar offset in this formu- lation, ensuring correct deviatoric stresses is sufficient to satisfy the momentum equation. Consequently, our training objective can be written:
[())
+ LBC + LIC Ldatas
where LBC and LIc enforce boundary and initial conditions, and Ldata penalizes any available labeled measurements. All loss terms are formulated in the least squares sense.
Cylinder Wake We simulate two-dimensional steady vortex-shedding flow at Re 100 around cir- cular cylinder centered at (0, 0). The computational do- main is defined as [L. 8] 2- 2] with the time interval [0, 7]. discretized in increments of At = 0.1. The do- main discretized using 40,000 interior collocation points to enforce the PDE residuals and 5,000 bound- ary and initial condition points apply the necessary constraints. A 4-layer, 50-neuron MLP with Tanh ac- tivation functions is trained using 50,000 iterations of the L-BFGS optimizer without any labeled data. Validation is performed using Direct Numerical Simu- lation (DNS) data from Raissi et al. (2019a) to assess the model's ''s performance. Comprehensive setup details and error analyses presented Appendix B.3.
Experimental Setups. We validate our approach on three representative incompressible Navicr-Stokes scenarios:
3D Beltrami Flow We consider the dinee-iimnsiooal Beltrami flow Reynolds number of Re to verify the accuracy of our RTNN framework. The computational domain is discretized using 2,601 inte- rior collocation points to enforce the PDE residuals, supplemented by 961 boundary and initial condi- tion points to impose the necessary constraints. An MI.P with hidden layers and neurons per layer employing Tanh activation functions is utilized. The model is trained using 100,000 iterations of the L- BFGS optimizer without any labeled data. Validation is performed against 26,000 interior points sampled within the domain to assess the model's performance. Detailed setup information and error plots are provided in Appendix B.2.
Results and Discussion. Table 3.3 summarize the perfor- mance of our RTNN formulation compared to NCL and standard PINNs by showing median average relative 7,2 er- ror across all fields and wall times, while Figures 4 and 6 summarize the training dynamics for these experiments and Figure 5 provides a qualitative comparison for the cylinder wake. As observed the table, our ansatz delivers improved
Steady Flow around a NACA Airfoil This experi- ment investigates the steady laminar flow at Re 1000
7

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Cyllnder
Airfoll
Beltrami
pressibility and momentum conservation in the DSFT Form:
Method
r1.2 Error
Time 65)
r1.2 Error
Time (s)
Error
Time (6)
RTNN
5.700c403
1,44-012
1.10e403
4.28c-44
297c402
NCL
2.54c-02
2-46c403
1.53-01
2,39c-43
1.73e-05
1.00c403
u 8 u T)
(21)
PINN
2.96442
3.12e402
2.48c.01
1.06e403
1.1640
1.82c402
Table 2. Comparison of RTNN, NCL. and PINN across three In- compressible Navier-Stokes experiments (Cylinder, Airfoil, Bel- trami). We report median relative 1,2 (rl.2) error and median wall times accross different seeds.
where u E R2 is:
T pl + u (Vu + (Vu))) + Madr
with Mader the the Maxwell magnetic stress.
Maev
B 2 I
B 6 B),
accuracy while achieving comparable training times. Across all three test cases, the RTNN approach consistently yields lower relative errors, indicating its potential for robust, data- efficient modeling of incompressible Havier-Stokes flows for both self-supervised learning in the PINN manner and in in scarce-data scenarios. Furthermore, for the cylinder wake, Figure 5 visually demonstrates RTNN's superior ca- pability in capturing the intricate vortex shedding structures. While the relative I error provides global measure of accuracy, the qualitative results in the figure highlight that RTNN yields solutions that are particularly effective in the complex wake region, where capturing fine details is crucial, an aspect that can be averaged out in global error metrics.
for magnetic field B.
RTNN Parametrization for 2D Incompressible MHD Let Sg be a (2+1) x (2+1) RTNN. From its block structure, we extract the physical fields as:
(ua, vo) (So)1:2.0
We parametrize the magnetic field B, via a separate network outputting scalar potential thg:
Bg = V (0, 0, vo) = (0y 0 R0).
V-B, = 0.
Define the Maxwell stress (including its isotropic part):
Ma = | - -Bo @ Bo).
3.4. Magnetohydrodynamiic (MHD)
We then let
We next consider the incompressible resistive magnetohy- drodynamics (MHD) equations, which couple fluid velocity and pressure magnetic field. The governing PDEs on domain Ω x [0, T] are:
og (S0) (uo, vp) @ (ug, vo) Ma.
ag is the purely fluid portion of the stress advection and magnetic terms have been subtracted. We then proceed in the same manner as section 3.3.
0,u + (u. v)u= -
( p
BI + +(B VJB +
Suu
σg = Pe 1 + oey, pa 1 ttr(os).
(18)
and define the viscous stress term
0,B + (u. 7) B B= -B=(B. V)u BB
(19)
Do = v (Vuo + (Vuo)
V =u=0,
BB
(20)
Training Objective. Enforcing momentum balance then requires ofev xe Deler We form residual
where (u. v) is the velocity field (2D case), p the
pressure, B = (Bz, By) is the magnetic-field vector, v is
the kinematic viscosity, is the magnetic diffusivity, and
Ra,visc 00
dev
Da-
2 = B + B2.
penalized in the loss. In addition, we include the induction equation residual
Piitto-Poteetial Formulation of the Magnetic Field. To enforce V. B- = 0 exactly, we parametrize B via a vector potential A:
Rinduction 3,Bg + (ug V) Bg (By () ug 17 V
Leading the final training objective:
B x A.
[(() R visc 122,
||Rinduction 112, LBo + LIc
In 2D, we may simply take A (0. 0. (xx
xr y)), that B = (0,th, - -.,b) automatically satisfies V B = 0. The induction equation (19) then becomes an evolution for vh:
Ldata
The first two terms penalize violations of momentum con- servation (viscous + magnetic) and induction equations, respectively. LBC and Lic enforce boundary and initial conditions, and Ldata integrates available labeled observa- tions. All loss terms are formulated in the least squares sense.
b.b + (u (..-aa tb = b.
Divergence-Free Symmetric Tensor (DFST) for Momen- tum. Similar the Navier- Stokes case, we unify incom-
8

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Reference Solution
Error (RTNN)
Error (Curl-Spinn)
0.125
0.8
0.100
0.075
0.4
0.4
0.050
0.2
0.2
0.000
0.25
0.50
0.75
x
Figure 3. Qualitative comparison of predicted velocity fields for the 2D periodic incompressihle MHD experiment (Re 1000, ReM 1000) at time t = 0.55. Left: Reference solution showing the velocity magnitude ||/w|||. Middle: Absolute error of the RTNN prediction, calculated as e ||uRed uRTNN Right: Absolute error the SPINN baseline prediction, calculated as URee USPINN ||2. This visualization highlights RTNN's 's enhanced accuracy in capturing the complex magneto- fendc dynamic structures compared to the SPINN baseline at the specified time.
Method
FL Error (Velocity)
rL, Error (B)
Wall Time (s)
qualitative comparison of the predicted velocity fields is pre- sented in Figure 3, which visually underscores the enhanced accuracy of RTNN in capturing the complex magneto-fluid dynamic structures compared the Curl-SPINN baseline. Our RTNN significantly improves velocity accuracy com- pared to the baselines and also enhances the accuracy of magnetic field predictions. Additionally, we demonstrated that RTNN can be incorporated into coupled systems for more complex problems.
RTNN
2.34e-02
1.04e-01
2340.41
CurI-SPINN
1.61e-01
1.55e-01
19
SPINN
1.8c-01
2.93c-01
414.73
Table 3. Comparison of methods for the MHD experiment, report- ing the median relative L2 error for velocity and magnetic fields in the table), and median wall time across five independent training runs with different random seeds.
Experimental setup We train and validate our method on a three-dimensional periodic incompressible MHD flow in j0, 2 with Reynolds numbers Re = 1000 and ReM = 1000. The simulation covers the time interval E 0.5] for training and tests at t = 0.55. Training data is generated using a spectral solver, with initial velocity un and magnetic field Bo sampled from Gaussian random fields. Periodic boundary conditions strictly enforced in all cases using the approach described in (Dong & Ni, 2021). We employ Separable Physics-Informne
4. Discussion and Conclusion
We introduced Riemann Tensor Neural Networks (RTNNs), new class of neural architectures tailored for encoding divergence-free symmetric tensors (DFSTs). By construc- tion, RTNNS exactly satisfy the DSFT conditions that en- codes conservation of mass and momentum. Our theoretical results confirm that RTNNs are universal approximators of DFSTs, and numerical benchmarks illustrate that they consistently improve accuracy compared to baselines such as standard PINNs and methods enforcing only part of the conservation (e.g., mass alone).
Neural Network (SPINN) (Cho et al., 2023) comprising 5 hidden layers and 500 neurons per layer to handle the structured 3D grid, discretized into 101 x 128 x 128 points. The model is trained using 50,000 iterations of the L-BFGS optimizer. Validation at 0.55 utilizes spectral solver data to evaluate performance. We compare three approaches: (i) our RTNN-based method, (ii) SPINN baseline with penalized residuals, and (iii) a Curl-SPINN parametrizing velocity and magnetic fields as the curl of a scalar potential. Additional details and performance plots are provided in Appendix B.5.
Limitations and Future Work While our experiments focus on fluid dynamics applications using MLPs, RTNNs have the potential to extend to diverse systems like Eu- Feu-ieu relativistic Euler, Boltzmann, to name but few. Future work includes enhancing algorithm performance, extending RTNNs to operator learning frameworks, and re- ducing computational costs by developing architectures with more efficient differential operators. Additionally, exploring function space optimization techniques (Jnini et al., 2024b; Jnini & Vella, 2025) could further improve the accuracy of RTNNs.
Results and Discussion Table 3.4 reports the relative L errors for both velocity and magnetic fields and Figure? shows the error evolution through training accross seeds. A

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Acknowledgments
Greydanus, S., Dzamba, M., and Yosinski, J. Hamiltonian neural networks. Advances in Neural Information Pro- cessing Systems, 32, 2019a.
A.J. acknowledges support from a fellowship provided by Leonardo S.A.A. This work was partially funded under the NRRP, Mission Component Investment 1.4, by the Euro- pean Union NextGenerationEU (proj. nr. CN 00000013).
Greydanus, S., Dzamba, M., and Yosinski, J. Hamiltonian neural networks, 2019b. URtpttps: https///axivvvoo abs 19699996.... 01563.
Impact Statement
Griffiths, D. J. Introduction to electrodynamics. Cambridge
University Press, 2023.
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which feel must be specifically highlighted here.
Haghighat, E., Raissi, M., Moure, A., Gomez, H., and
Juanes, R. A physies-informed deep learning framework for inversion and surrogate modeling in solid mechanics. Computer Methods in Applied Mechanics and Engineer- ing, 379: 1 13741, 2021.
References
Bonfanti, A., Bruno, G., and Cipriani, C. The challenges of the nonlincar regime for physics-informed neural net- works. In The Thirty-eighth Annual Conference on Neu- ral Information Processing Systems, 2024, URL https: //openreview wtttleeweeeeeeeeeeeeeeeee wwe//fftroooeeeeeeee forum? d== 6vPt.TT....
Hendriks, J., Jidling, C., Wills, A., and Schön, T. Lin- early constrained neural networks. arXiv preprint arXi::200..0160) 2020.
Hu, Z., Shukla, K., Karniadakis, G. E., and Kawaguchi, Tackling the curse of dimensionality with neural networks. arXiv preprint urXiv:2307. 12306, 2023.
K.
Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D.. Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q. JAX: composable transformations of Python+NumPy programs, 2018. URL http::/ i thub. c googl 1 e jax.
Jnini, A. and Vella, F. Dual Natural Gradient Descent for Informed Neural Networks.
Scalable Training of Physics-Iffrrmed arXiv preprint arXiv:2505. 1404, 2025.
Cai, S., Mao, Z.,, Wang, Z.,, Yin, M., and Karniadakis, G. E. Physics-informed neural networks (PINNs) for fluid me- chanics: A review. Acta Mechanica Sinica, 37(1 1738, 2021.
Jnini, A., Goordoyal, H., Dave, S., Korobenko, A., Vella, F,, and Fraser, K. Physics-constrained deepONet for surro- gate CFD models: a curved beckward-facing step case. ICLR 2024 Workshop on A14DifeerentallEquation In Science, 2024a. URL https: / opwnreveeeeeeeeeeeeeeeeee .nt//
Cho, J., Nam, S., Yang, H., Yun, S.-B., Hong, Y. and Park, E. Separable physics-informed neural networks, 2023.
URL /ettps://arxtae org/ /bs/abba/ 11969,
Jnini, A., Vella, F., and Zeinhofer, M. Gauss-Newton Nat- ural Gradient Descent for Physies-Informet Computa- tional Fluid Dynamics. arXiv preprint arXiv:2402. 2024b.
Cranmer, M., Greydanus, S., Hoyer. S., Battaglia, P., Spergel, D., and Ho. S. Lagrangian neural networks, 2020.
URL https: arxiv. org/ abs 2003. 04630.
Karniadakis, G. E., Kevrekidis, I. G.,, Lu, L., Perdikaris, P., Wang, S., and Yang, L. Physics-informed machine learning. Nature Reviews Physics, 3(6):422-440, 2021.
Dong, S. and Ni, N. A method for representing periodic functions and enforcing exactly periodic boundary con- ditions with deep neural networks. Journal of Compu- tational Physics, 435:1 June 2021. ISSN 0021- 9991. doi: 202. 10244. URL http:// dx www.oo....oooo.... org 10. 1006/1 jep 2021. 110242.
Kidger, P. and Lyons, T. Universal approximation with deep narrow networks, 2020. URL htts:t/:://rxx//oorr
abs/1905.0855 08539.
Ethier, C. R. and Steinman, D. A. Exact fully 3d solutions for benchmarking. International Journal for Numerical Methods Fluids, 19:369-375, 1994. URL https: htt//app..aaaaaaee http://///appssmmattcccooooo/// sem.ssttcncccchl..rr org/Corpus
Kovachki, N., Li, Z., Liu, B., Azizzadeneshelli K., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Neural operator: Learning maps between function spaces. arXiv 2021.
preprint arXiv:2108.
62789476.
Lagaris, I. E.,, Likas, A., and Fotiadis, D. I. Artificial neu- ral networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 9(5): 987-1000, 1998.
Giles, C. L.. and Maxwell, T. Learning, invariance, and generalization in high-order neural networks. Applied optics, 26(23):4972-499 4978, 1987.
10

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. Gradient- based learning applied to document recognition. Proceed- ings of the IEEE, 86(11):228 86(11):2278-2324, 1998.
Wang, S., Teng, Y., and Perdikaris, P. Understanding and mitigating gradient flow pathologies in physics-infommn neural networks. SIAM Journal on Scientific Computing, 13(5)::3055A3300 2021.
Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhat- tacharya, K., Stuart, A., and Anandkumar, A. Fourier neural operator for parametric partial differential equa- tions. XrXiv preprint arXiv:2010.00 08895, 2020.
Liu, N., Fan, Y., Zeng, X., Klöwer, M., Zhang, L.,, and Yu, Y. Harnessing the power of neural operators with automatically encoded conservation laws, 2024. URL atttps://arxiee
abs abs/2312.
Müller, J. and Zeinhofer, M. Position: Optimization in sciml should employ the function space geometry, 2024. URL https: aarxivv.. org/ www./44/4000 07318.
Nocedal, J. and Wright, S. J. Numerical optimization. Springer, 1999.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physics- informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019a.
Raissi, M., Perdikaris, P., and Karniadakis, G. E. Physics- informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-7017 2019b.
Richter-Powell, Lipman, and Chen, R. Neural conservation laws: A divergence-free free perspective. Ad- vances in Neural Information Processing Systems, 35: 38075-38088, 2022.
Serre, D. Solutions classiques globales des équations d'euler pour un fluide parfait incompressible. Annales de I'Institut Fourier, 47(1):139-153. 1997. URL
https: /ww.. nwwdmm htt.mddm..og//rticlessssss
5802/a1f.1565
Serre, D. Divergence-fin positive symmetric tensors and fluid dynamics. Annales de l'Institut Henri Poincaré 2018. URL (Analyse Non Linéaire), 1705. 00331.
https://arxiv.. wttp//.rrg///aa. abss
Serre, D. Compensated integrability. applications the vlasov-poiee equation and other models of mathemat- ical physics. Journal de Mathématiques Pures Ap- pliquées, 127:67-88. 88, 2019.
Serre, D. Hard spheres dynamics: Weak vs strong collisions. Archive for Rational Mechanics and Analysis, 240(1): 243-264, 2021.

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
A. Proofs
A.I. Proof of Theorem 2.1
In this section, we prove the main result on the representation of divergence-free symmetric tensors (Theorem A.I), which states that ab (symmetric and divergence- free)) exactly coincides with the image of the map Kacbd H> V "VV"Kachdl where Kaebd satísfies Riemann-like symmetries.
We first recall the classification of Riemann-like (0, 4)-tensors (Lemmaa A.2), and then prove the surjectivity of the map Φ: Kacha + V Kaebd (Lemma A.4). Combining these two ingredients completes the proof of Theorem A.I.
A.1.I. MAIN THEOREM: REPRESENTATION OF DIVERGENCE- FREE SYMMETRIC TENSORS
Theorem A.1 (Representation of Divergence- Symmetric Tensors on a Flat Manifold), Let be an n-dimensional real vector space with a fixed basis {ea (e and let {er n denote the corresponding dual basis of V*. Let A² denote the space of 2-forms on V. Consider the space of all (0, 4)-tensors Kabed defined on a flat manifold equipped with a Levi-Civita connection V, satisfying the following symmetries:
1. Anfisymmetry within index pairs:
(aabyod 0, Kab(od) 0,
2. Symmetry between pairs:
Kabod Kedab.
Let (wi....,wm) be a fixed basis of A2V*, where 2
m(n-1),
Define the tensors:
T(is) wi(t A wnle A Aez) wy(e A ef) e2)
Then:
1. The space of divergence-free, symmetric (0, 21tensoors Sab on the flat manifold is exactly the image of the map
Kachd b Sab =
V Kaebd-
2. Moreover, any such Sab can be expressed as
Sab
This)) acta
Cijs
Isisism
where Cij are smooth scalar functions.
Proof Outline. By Lemma A.2, every Kachd with the above symmetries can be expanded in the basis {T(i.)). Hence 7"""KKasa can be written in the claimed form. We also check that VavaKantt is symmetric and divergence-fir using flatness (commuting covariant derivatives) plus the antisymmetrics. Finally, Lemma A.4 shows that any symmetric, divergence-fres Sab arises from some Kaetd, establishing that the image of the map is the space of such Sab-
A.1.2. CLASSIFICATION OF RIEMANN-LIKE TENSORS
Lemma A.2 (Classification of Riemann-like (0, -ennors) Let V be n-dimensional real vector space, and let A2V* be the space of 2-forms on V. Denote
m = dim(A?V*) nen - 1)
2
Consider the vector space
T Tabed
Tabed
-haad:
Tabed -T Tabdez
Tabedt Tedab) )
12

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
There is canonical veclor-space isomorphism between this space and Sym (A2V*). In particular, its dimension is
n(n-1) n(n-1) 1)
+1)
m (m +
2
2
2
Moreover, if is basis for V then corresponding basis in the space of such (0, -eesssrs is given by
The
(w; Wy wy + wi @ wi) (ez A e nes), isisism
abidd
Proof of Lemma A.2. Step 1: From T to a symmetric bilinear form. Given Tabed with the stated symmetries, define for 0,3 E A²V*:
a ab ea A c 3 3od e A Edi
T(a. 3) = Tabod ab god.
Since ab and god are antisymmetric in their respective index pairs, and Tabed is antisymmetric in (a, b) and (c, d), this well-defined. The symmetry Tabest Todab implies that T symmetric as bilinear form on A²V*. Hence T Sym Sym?(A3V") (A3V*).
Step From a symmetric bilinear form back to Conversely, given T Sym?() define
Tabed 7(% A cbs e A eÿ)
A straightforward check shows that Tabed inherits antisymmetry in (a, b), antisymmetry in (c, d), and the pair-exchange symmetry (a, b) ++ (c, d).
Step 3: Isomorphism and basis. These two constructions are linear inverses of each other, yielding a vector-space isomorphism between our (0, 4)-tsssosss and Sym Sym2 (A?V*.. The dimension follows from standard linear algebra. (.i,.... 4mm } is basis of A2V*, then {wi & wj + wj @ wi : Kiisj S m) forms basis for SA3V() Mapping these to (0, 4)-tensors via the above correspondence yields the stated basis (T(W))
A. 1.3. SURJECTIVITY OF THE MAP Φ: K A.AKK
As discussed, the second key ingredient for Theorem A.1 is showing that every divergence-free, symmetric (0, 2)-tensor Sab can be obtained from some Riemann-like (0, 4)-tensor Kuedd. Equivalently, the map
+ Kacbe
Sab
V Kaebd
is surjective.
Lemma A.3 (Surjectivity of 4). Let Kaeba be any Riemann-like (0, 4)-tensor (i.e. satisfying the symmetries of Lemma A.2). Define S.b Karhtt Then Sab is automatically symmetric and divergence-free. Moreover, 4 is onto: for any given symmetric, fivergence-free (0, 2)-tensor Sob. there exists Riemann-like Kachd such that Sab = V K
Kaeha.
Proof. Part 1: We first verify that if Sab = - Kaebe, then Sab is symmetric and divergence- free.
SSmmmery: Using Kaetd = Kidac, we get
Sab = A'Vd Kaete D Kodac = Sba-
Diverrencc--fre:: flat space, covariant derivatives commute, so
V Sa. = 4'A. K acbil d Kactd-
Because Kacha is antisymmetric in (a, e), one shows that V K Kaetd = Hence S 0.
Part 2 (Surjectivity):: We now show that given any symmetric, divergence-frre Sas, we can solve Sab = V'V Kachd for some Riemann-lit like Kacdd-
13

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
A.2. Surjectivity of the map Theorem A.4 (Surjectivity of the Q). We will now prove that the following map surjective:
k S
(22)
Sab = $ (Kahd)) = VO K
(2.)
where Kacha is Riemann-like (0, 4) tensor and Sab is a symmetric divergence-free (0, 2) tensor.
Proof. This is equivalent to showing that for any given Sab. there exists Kaeha such that San - 4 (aaaohd). As in previous sections, we assume we are working on flat manifold.
Let's define Poisson's 's equation (Griffiths, 2023) for each index of the tensor Sab. with appropriate boundary conditions desired:
Laa Sab
(24)
These Poisson's equations can be solved by the well-known method of Green's function (Griffiths, 2023).
By solving the previous equations, we obtain a tensor Lab with the following properties:
Lab = Lba
(25a)
4 L.ab
(25b)
Sab ALab
(25c)
It's easy to show that:
Eq. 25a follows from the fact that the Laplacian operator A in Eq. 24 is meant componentwise and Sab is symmetric.
Eq. 25b follows from the fact that the derivatives commute on a flat manifold.
Eq. 25c follows form definition 24.
We now easily define tensor Kaeba with the desired properties from Lub- Just let:
Kaetat 2 [S Ldja)
(26)
= Sab SaaTsc TetT + Sed Taba
We now verify the Riemann-like properties of K:
K (accba = 4 (ōapILaf - Sopt + 8op Lajo 5
(27)
Kac(ba) 4 (dajoLetcc SeppLaaa + dajaLijc dejalibja)
(28)
Kacha 2 (dajoLedec- Laja)
(29)
babLide duaLic dibadda + Seabiba
Sea Led dao Las Lad + d Lab
Sma Lood - 6nc Load dea Doo + 8 Taab
2 (ōota Ioja 8a Lage)
Kidac
VO K acda
2 A'Aa (dusoLafc Sc)bLida))
(30)
A'A Lan) V-Va (5ag D-Ad (b.L + A-AA (Soabma)
San V'Va - V. Line deb V. + Saa V'rdd
- 0- 0 V'VeLba
ALab
Sab
14

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrained Networks
So, since we build inverse for every Sab. the map Φ Kaebd + Sub is surjective.
A.2.1. COMBINING THE LEMMAS TO PROVE THEOREM A.1
Proof of Theorem A.I. (1) Symmetry and divergence-free for V Kacbd. By Lemma A.4 (Part 1), the tensor yopa aebd is always symmetric and divergence-free if Kachd is Riemann-like.
(2) Surjectivity: Every symmetric, divergence-free Sab arises from some Kaohd- By Lemma A.4 (Part 2), the map
@ + V'V is onto the space such Sab.
(3) Basis and explicit representation. From Lemma A.2, any Riemann-like Kacha can be expanded the basis (T())). Consequently,
V.V WV7dKabbd
Kaebd
ThE acha
V'V
1515i5m
for some scalar functions {aj). This shows that every Sab can indeed be written in the claimed form.
Hence all parts of the statement hold, and the proof is complete.
A.3. Proof of theorem A.5
Theorem A.5 (Universal Approximation for RTNN), Ler Ω C R be compact domain, and let S E C(2, R3x3) be any divergence-rree symmetric tensor field on Q. For any > O, there exists a Riemann Tensor Neural Network So with a fixed narrow width and arbitrary depth, such that
sup |[((x) - Sx(x)) <G
(31)
where || denotes the Frobenius norm.
Proof. To prove Theorem A.5, we proceed two main steps:
1. Surjectivity of the Map from K to S: By Theorem A.4, the map
Φ KK + S,
(32)
defined by
Sab = V Kaehds
(33)
is surjective. This implies that for any divergence-free symmetric tensor field S E C(Ω; Rmxm), there exists a tensor with Riemann-like symmetries such that
5=
(34)
2. Approximation Tensor K Using Neural Networks: Since K a linear combination of scalar functions, specifically
Kaebd Σ (x)Tr))))
(35)
where T((..) are fixed basis tensors and ciy (x) are scalar coefficient functions, we can approximate each cty (z((x) using deep narrow neural networks.
By the Universal Approximation Theorem for deep narrow neural networks (Kidger & Lyons, 2020), for each scalar function (z) and for any >> 0, there exists neural network NNo., such that
sup (x) =
city (x;; 0zj)
(36)
15

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Plysies-Constrainee Networks
where C' is constant dependent the norms of the basis tensors T(5) and the chosen tensor norm || Construct the approximated tensor K as
||.
Kx(x)
Σ 2
NN cj (x; 0.y (Ti)))
(37)
Applying the map & to Kg, we obtain the approximated tensor So:
So(z) (<))(x)) Ko(z).
(38)
Error Propagation:
The approximation error in each cy (x) propagates linearly to (a))). Specifically, we have
[5(x) - So(x)]
(39)
V'V - [K])
A44a
Σ (cy(()) (a;
(x; (zg)) 1) actid
|
sup A'V" (G)(x) - ci (x; 0z))) BE
|لااااي||
sup (ي)()) 2 (z; 01s) Cis
| C
where Cg accounts for the bounds on the second-order partial derivatives V'V(ay(z) over the compact domain Ω and the final equality holds by the choice of C.
Therefore, by appropriately choosing the neural networks NNe., to approximate each scalar coefficient (1y(2) within e/C, we ensure that the approximated tensor (1(x) satisfies
sup |[(x)) - Sa(z)] <e senl
(40)
A.4. Equivalence of Explicit Momentum Residual and Isotropic Loss Term
Consider the compressible Euler momentum equation (neglecting viscous terms):
(.(p)) - V (pu @ u + pl) = 0,
where p is the density, u the velocity field, and p is the pressure.
In our RTNN parametrization, the network outputs flux tensor defined as:
( p pu
(pu)T
S
pu
with p So,0 and pm S1,0- The stress tensor is then defined by:
= S1.2.1:2 p.u. 09 n.
16

---

Riemann Tensor Neural Networks:
Learning Conservative Systems with Physies-Constrainee Networks
We impose zero-deviatoric constraint on o, meaning its deviatoric part vanishes:
1 -tr(()) ( I= =0
dev
This forces:
1 tr(a) 1 pl,
= n
With the pressure defined as:
p tr(o).
Substituting p7 back into the flux tensor, get:
(pu)T S pu pu E u pT)
Our architecture enforces that S is divergence-free:
(a)(m)) V (pu)
S = (pm) V (pu pT)
The first row yields the continuity equation, which vanishes:
Oept V (pu) O.
The second row yields:
0.(pu) - (pu Q u = 0,
which, when expanded, yield the following equations:
0.p + V (pu, pv) = 0
(41)
Op(p)) - V (pu², puv) -r.
(42)
(pv + (pue,po² -0P
(43)
This is exactly the system in (11)-113).
In other words, the condition on the stress tensor is numerically equivalent to having the momentum residual vanish. A similar reasoning can be applied when adding viscous terms (or magnetic-terms)) for instance.
Furthermore, to prove our mathematical derivation numerically, we conducted an experiment in which the velocity and pressure fields extracted from single RTNN were trained using the classical momentum equation rather than the simplified condition the stress tensor. In this setting, the only difference among the experiments is the underlying architecture. Although this process is more computationally expensive, owing to the necessity of computing a fourth-order derivative, the RTNN-trained fields still demonstrate similarly superior accuracy compared to those obtained with PINN and NCL. We summarize the obtained results below for the two unsteady Navier-Stokes problems: flow around a cylinder and the Beltrami flow (see Table 4).
PINN and NCL were both trained using the full system (11)-(14), which includes the conservation of mass, momentum, and energy equations. For NCL, the mass conservation equation is satisfied by construction and cancels up to machine precision; similarly, for RTNN, conservation of mass is similarly enforeed structurally and cancels out. The energy equation (14) is also included in the loss, similarly to the other methods.
The only apparent difference is that RTNN replaces the explicit momentum residual with an isotropic loss term. However, this isotropic loss term naturally emerges when substituting the ETNNN-predict fields into the flux-form momentum equation. Imposing the zero-deviatoric condition on the stress tensor is numerically equivalent to enforcing the momentum equation directly in its expanded form.
17

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Beltrami
Cylinder
Method
L2 Median
L2 IQR
L2 Median
L2 IQR
RTNN (extended momentum equation)
4.55e-04
8.1e-05
5.13e-03
4.2e-04
RTNN (stress tensor)
4.28e-04
5.9e-05
5.70e-03
4.3e-04
NCL
1.74c-03
8.5c-05
2.54c-02
1.9c-03
PINN
1.41e-03
3.1e-04
2.99e-02
7.9e-04
Table 4. Comparison of methods across Beltrami and Cylinder tasks. Reported: median and IQR of relative L2 error.
B. Additional resources for Experiments
B.1. Euler Equation
The simulation of the 2D isentropic Fuler vortex is based analytical solutions, ensuring precise modeling of vortex dynamics. Below, provide detailed initial and boundary conditions, well summary of the experimental setup parameters and hyperparameters.
Table 5. Setup and Hyperparameters for the Fuler Vortex Experiment
Parameter
Value
Optimizers
BFGS
Architecture
4-layer Multilayer Perceptron, width 50
Activation Function
Tanh
Domain
Ω [0, [.) x [0, Ly]
Time Interval
[0,
Collocation Points (Interior)
1,000
Collocation Points (Boundary)
200
Validation Points
10,000
Evaluation Metric
Relative _ Error
The initial density, velocity components, and temperature distributions are defined based on the analytical solution of the Euler vortex:
Initial Conditions:
T(x,y)
1/(7-1)
p(z, 0, = Poo
(44)
Too
u(x, y, 0) u
Uc) exp [1 ,² v)].
(45)
2T
o(x,y,00 0) = Uoa +
(xr Ie)) exp [1 , v)
(46)
2x
(7 1)
T(x v) Tox
8x exp [2(1 7 w))]
(47)
where 72 y) = - (x)= (y y. represents the radial distance from the vortex center (xe: Be).
Boundary Conditions: On the boundaries of the spatial domain 00, we impose fixed values for density, velocity, and pressure:
plan = P3os (u, v)aa = (Moc, Voo), plan KP
18

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
B.2. Beltrami Flow
Table 6. Setup and Hyperparameters for the Beltrami Flow Experiment
Parameter
Value
Optimizers
BFGS
Architecture
4-yayer Multilayer Perceptron, width 50
Activation Function
Tanh
Domain
Ω
1 x 1, [-1,1
Time Interval
[0.1]
Collocation Points (Interior)
10,000
Collocation Points (Boundary)
961 per face 10,000
Validation Points
Evaluation Metric
Relative 1,2 Error
We consider the unsteady three-dimensional Beltrami flow originally described by Ethier & Steinman (1994) at Re 1. The spatial domain is Ω = [-1, 11s and time ranges over [0, 1]. The exact solutions are:
u(x,8,,,1, = = sin(y + 2) + e² cos(æ + y) -
o(2,8,211) = -e sin(z + x) + et cos(y + 2)
02,, = -e" -insin(x + y) + et cos(2 2 + z) e-
and
= - - + 22 + = + 2sin(x + y) cos(2 + x)) c
+ 2yinny + +) cos(x + y) -+ + 2sin(2 + x) cos(y + 2) c) 2-1
These expressions satisfy the incompressible Navier-Stokes equations exactly, making the Beltrami flow an ideal test for
verifying numerical methods. We use 10,000 interior collocation points to enforce the PDE residual and 961 boundary/initial points per face for constraints; an additional 10,000 points in the interior serve as a validation set. We report the training dynamics in Figure B.2.
19

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
RTNN
NCL.
PINN
109
Relative L2 Error
10
10-2--
10-1
0
20
40
60
80
100
Simulation Time (s)
Figure Comparison training relative 2 error over simulation time for the 3d Beltrami experiment across methods. The solid lines represent the median relative L error, while the shaded regions indicate the interquartile range (IQR). RTNN achieves the lowest L2 error, followed by PINN and NCL. The x-axis represents the simulation time seconds), while the y-axis shows the relative L error on logarithmic scale.
B.3. Cylinder Wake
Table 7. Setup and Hyperparameters for the Cylinder Wake Experiment
Parameter
Value
Optimizers
BFGS
Architecture
4-layer Multilayer Perceptron, width 50
Activation Function
Tanh
Domain
Ω [1. x 2- 2]
Time Interval
[0. 7], A: = 0.1
Interior Collocation Points
40,000
Boundary/Initial Points
5,000
Evaluation Metric
Relative L Error
For the two-dimensional cylinder wake at Re = 100, we consider circular cylinder centered at 0). The downstream domain is Ω = [1, 8] x [-2, 2] in space, with the simulation time interval [0, 7| discretized in increments of At = 0.1. We train a 4-yeyrr MLP (width 50, Tanh activations) on 40,000 interior collocation points and 5,000 boundary/initial points, without using labeled data. We validate the solution against direct numerical simulation (DNS) from Raissi et al. (2019a). We report the training dynamics Figure 5
20

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
RTNN
NCI.
PINN
109
Relative L2 Error
10-1
10-2
0
10000
20000 30000
40000 50000
Iterations
Figure Comparison training relative L error over simulation time for the Cylinder wake experiment across all methods. The solid lines represent the median relative L error, while the shaded regions indicate the interquartile range (IQR). RTNN achieves the lowest L2 error, followed by NCL and pinn. The x-axis represents the number iterations, while the y-axis shows the relative L error on logarithmic scale.
B.4. NACA Airfoil
Table 8. Setup and Hyperparameters for the NACA Airfoil Experiment
Parameter
Value
Optimizers
BFGS
Architecture
4-layer Multilayer Perceptron, width 50
Activation Function
Tanh
Evaluation Metric
Relative 2 Error
We investigate the steady laminar flow at Re = 1000 around a NACA 0012 airfoil. This problem is recast as steady Navier-Stoker system by treating time as a dummy dimension. An MLP with 4 hidden layers (width 50, Tanh activations) is trained with labeled data from our in-house RANS solver with 2'000 data points and 40'000 collocation points for PDE residuals. We use 12'000 data points to validate the recovered fields. We report the training dynamics in Figure 6
B.5. MHD
train and validate our method on a three-dimensional periodic incompressible Magnetohydrodynamics (MHD) flow within the domain [0, 2)3 with Reynolds numbers Re 1000 and ReM 1000. The simulation covers the time interval t 10, 0.5] for training and initial velocity un and magnetic field By sampled from Gaussian random fields. Periodic boundary conditions strictly enforced to maintain physical consistency across the domain. We employ Separable Mayssics-Iffrmmen Neurul Network (SPINN) (Cho et al., 2023) with hidden layers and 500 neurons per layer to handle the structured 3D grid, discretized into 101 x 128 x 128 points along each respective dimension. The model is trained using 50,000 iterations of the L-BFGS optimizer to effectively minimize the loss function. Validation at = 0.55 utilizes spectral solver data to assess performance and ensure accurate replication of MHD flow characteristics. We compare three approaches: (i) our -TNN-bassd method,
0.55. Training data is generated using a spectral solver, ensuring high accuracy, with
21

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
RTNN
NCL.
PINN
10a
Relative L2 Error
10-1
0
20000
40000
60000
80000
100000
Iteration
Figure Comparison of training relative error over simulation time for the NACA Airfoil experiment across all methods. The solid lines represent the median relative L error, while the shaded regions indicate the interquartile range (IQR). RINN achieves the lowest L2 error, followed by PINN and NCL. The x-axis represents the simulation time (in seconds). while the y-axis shows the relative L error on logarithmic scale.
(ii) SPINN baseline with penalized residuals, and (iii) a Curl-SPINN parametrizing velocity and magnetic fields as the curl of scalar potential, Additional details and performance plots are provided below.
Table 9. Setup and Hyperparameters the MHD Experiment
Parameter
Value
Optimizer
L-BEGS
Architecture
Separable Physics-Informed Neural Network (SPINN)
Hidden Layers
5
Neurons per Layer
500
Activation Function
Tanh
Domain
[0, L]3
Time Interval
Training: & E 0.5), Testing: / 0.55
Grid Discretization
101 128 x 128 points
Reynolds Number
Re = 1000
Magnetic Reynolds Number
ReM 1000
Training Iterations
50,000
Validation Data
Spectral solver data 0.55
Evaluation Metric
Relative 2 Error
22

---

Riemann Tensor Neural Networks: Learning Conservative Systems with Physies-Constrainee Networks
Velocity Error
Magnetic Potential Error
RTNN
RTNN
CUIL
CURL
PINN
PINN
Relative L2 Error
10
20000
40000
60003
50000
100000
20000
10000
80000
100000
Iterations
Iterations
Figure Comparison training relative L error over iterations for the MHD experiment across methods. The left subplot represents the velocity relative L error, while the right subplot represents the magnetic potential relative L error. The solid lines indicate the median relative _ error, and shaded regions depict the interquartile range (IQR). RTNN achieves the lowest L error both velocity and magnetic potential, followed by PINN and NCL. The x-axis represents the iterations, while the y-axis shows the relative L error on logarithmic scale.
23


Reviews:
[{"review_id": "I8LbqaDsgw", "reviewer": "Reviewer_yLyj", "summary": "The paper proposes RTNNs, which can encodes the divergence-free constraints in neural networks within the PINN framework. The divergence-free constraint is satisfied by computing the hessian matrix of a feed-forward MLP network and combining with a special-designed basis. The method is evaluated on a range of PDE tasks requiring the divergence-free condition.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their time and feedback. Below we address the reviewer's questions:\n\n**Label Data Inclusion**  \n\nWhile we can incorporate labeled data to penalize discrepancies between predictions and observations, for this experiment we did not. We acknowledge the inconsistency between line 270 and line 236 and will correct it in the final version.\n\n**Both methods use similar MLP architectures for fairness. Are they exactly the same?**  \nBoth methods share identical hidden-layer configurations (depth, width, activation). The only difference is the *final output layer*, where each architecture produces the specific number of scalars required by its formulation. \n\n**Is NCL a Strong Baseline?**  \nWe address the NCL performance issue for the compressible case in our response to reviewer G2wL.\n\n**Will computing the Hessian matrix affect the scaling of computation cost?**  \nA key contribution of our work is the formulation of a computationally tractable ansatz that scales well in the 2D+1 and 3D+1 cases of interest. Rather than naively parameterizing a (0,4)-tensor— which would require outputting \\\\( \\mathcal{O}(N^4) \\\\) , where \nN is the number of spacial components components — we exploit an isomorphism with two-forms to reduce the number of scalar outputs to \\\\( \\tfrac{N(N+1)}{2} \\\\), thereby limiting the number of second-order derivatives that must be computed.\n\nComputing the Hessian is inherently expensive; however, we address this challenge through two derivative evaluation strategies:\n\n- **Taylor-mode Automatic Differentiation:**  \n  We compute second derivatives using Taylor-mode AD, which efficiently propagates truncated Taylor expansions (\"jets\"), thereby sharing repeated derivative computations across different orders through the network. For a scalar-output MLP, the complexity scales as \\\\( \\mathcal{O}(L d W^2) \\\\), where \\\\( L \\\\) is the number of layers, \\\\( W \\\\) is the layer width, and \\\\( d \\\\) is the order of the derivative. This ensures that, for typical 2D+1 and 3D+1 PDEs, the cost of higher-order derivatives remains within a modest constant factor of a standard forward pass.\n\n- **SPINNs for Magneto-Hydrodynamics:**  \n  In Section 3.4, we take advantage of the SPINN architecture[1] to decompose multi-dimensional derivative evaluations into per-axis components, reducing the cost from \\\\( \\mathcal{O}(N^d) \\\\) to \\\\( \\mathcal{O}(N\\ d) \\\\), where \\\\( N \\\\) is the number of collocation points per axis and \\\\( d \\\\) is the order of the derivative.\n\nTogether, these strategies ensure that while computing the Hessian is expensive, our approach remains scalable and computationally tractable for the 2D+time and 3D+time problems of interest.\n\n**What do the parentheses in equation 2 and 3 in the subscripts mean?**  \nIn tensor algebra, a pair of indices enclosed in round parentheses denotes symmetrization with respect to those indices. For example:\n\n\\\\[\nT_{(ab)\\dots} = \\frac{1}{2} \\left( T_{ab\\dots} + T_{ba\\dots} \\right)\n\\\\]\n\n**Intuition Behind Proof 2.1:**  \n\n At each grid point, we must generate a flux tensor that satisfies the DSFT condition exactly. Theorem 2.1 provides a systematic way to achieve this. We take advantage of the fact that the divergence of a second derivative of a four-tensor naturally vanishes under certain symmetries (2)--(4). Specifically, if we define $S_{ab}$ in terms of a second derivative of some four-tensor $K_{abcd}$, the divergence-free condition is automatically satisfied due to the symmetry properties of partial derivatives.\n\nThis is reminiscent of how potentials are used in physics to impose constraints on derived quantities. For example, in electromagnetism, introducing a four-potential \\\\( A_\\mu \\\\) ensures that the field strength tensor\n\n\\\\[\nF_{\\mu\\nu} = \\partial_\\mu A_\\nu - \\partial_\\nu A_\\mu\n\\\\]\n\nautomatically satisfies the Bianchi identity:\n\n\\\\[\n\\partial_\\lambda F_{\\mu\\nu} + \\partial_\\nu F_{\\lambda\\mu} + \\partial_\\mu F_{\\nu\\lambda} = 0\n\\\\]\n\nWe hope that we addressed the reviewer's questions sufficiently and we will stay attentive for further questions throughout the review process.\n\n**References :**\n\n\n[1] ​Cho, J., Nam, S., Yang, H., Yun, S.-B., Hong, Y., & Park, E. (2023). Separable Physics-Informed Neural Networks. Advances in Neural Information Processing Systems 36 (NeurIPS 2023).", "claims_and_evidence": "- Although understanding the theorem in section 2.1 seems to be beyond my knowledge, the description of divergence-free construction using hessian matrix is clear and convincing.\n\n- The experimental results look strong compared to PINN and NCL, although I am not familiar with the latter one.", "methods_and_evaluation_criteria": "- The benchmarked problems are reasonable and diverse.\n\n- In line 270 left it says \"$L_{data}$ incorporates supervised learning by penalizing the discrepancies between the model predictions and observed data labels\". But in line 236 right it says \"training is performed entirely without labeled data\". It seems to be inconsistent.", "theoretical_claims": "- I did not check the correctness of math but the theory seems to be sound.", "experimental_designs_or_analyses": "- Line 243 right: \"Both methods use similar MLP architectures for fairness.\" Are they exactly the same?\n\n- Is NCL supposed to be a strong baseline? The results of NCL in table 1 seems to be worse than PINN.", "relation_to_broader_scientific_literature": "- They are adequately discussed.", "essential_references_not_discussed": "- Not I am aware of.", "other_strengths_and_weaknesses": "- The paper is well written.\n\n- The method still requires other soft constraints to satisfy other constraints.\n\n- The use neural network model is rather simple.", "other_comments_or_suggestions": "- I did not notice typo.", "questions_for_authors": "- Will computing Hessian matrix affect the scaling of computation cost?\n\n- What does the parentheses in equation 2 and 3 in subscripts mean? \n\n- Is there any intuitive way to understand theorem 2.1? For example, what does it mean in terms of quantities on the simulation grid?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "OH0RT0ATWc", "reviewer": "Reviewer_pPMm", "summary": "This paper concerns preserving certain structure properties when using neural networks to solve partial differential equations. The method uses an inductive bias by enforcing a certain form of the tensor field approximation, and the final form is called \"Riemann Tensor Neural Network\" (RTNN) involves some non-trainable final layer, plus an \"differentiation layer\" where the Hessian of the output from the previous layers is computed.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the detailed comments and suggestions. We appreciate them and hope to address some of the questions below:\n\n**Suggestions:**  \nWe agree that including representative solution plots would help the reader infer the difficulty in PINNs solution. We also appreciate the suggestion to cite relevant literature on structure-preserving numerical methods. Both additions will be included in the camera-ready version.\n\n**Why do RTNNs perform better than PINNs:**  \n\nPINNs are notoriously difficult to optimize despite their recent success [1] . From the very advent of PINNs, several authors have advocated modifying machine learning models such that the boundary conditions or conservation laws are exactly zero. This has been empirically shown to improve training accuracy, including in our own work.  The work in **De Ryck et.al** [2] links the training dynamics of PINNs to the condition number of a Hermitian operator A that depend both on the PDE operator and the boundary operator. They show that hard-imposing constraints within the architecture reduces the condition number of A, which in turn improves the training dynamics and accelerates convergence.\n\nAnother way of explaining the superior performance of RTNNs relative to unconstrained networks is their ability to capture and exploit the intrinsic relationships between physical fields. In many fluid dynamics problems, one must simultaneously model a scalar density field, a vector velocity field, and a scalar pressure field. Conventional unconstrained models typically treat these components as independent scalar output channels, thereby ignoring the important geometric and physical correlations—such as the interdependence between the directional components of the velocity field, or the influence of pressure gradients on both density and velocity. By contrast, RTNNs treat these quantities as a unified field coupled through a shared conservation-law-based flux, this in turn leads to superior accuracy during the training. \nWe believe this discussion will improve the quality of our manuscript. We will include a dedicated paragraph on this topic in the final version, and we thank the reviewer for raising this point.\n\n**Question on Differentiation Layer:**  \nIn our implementation, we compute at most third-order derivatives, and all experiments were conducted in double precision (we will state this explicitly in the revised manuscript). While repeated differentiation can in principle lead to numerical issues, we did not observe any instabilities or degradation in constraint satisfaction. The resulting tensors remained divergence-free and symmetric to machine precision throughout. In fact, we also verified that this property holds under single precision as well (although naturally this leads to degradations in accuracy for all methods).\n\nWe hope that we addressed the reviewer's questions and we will stay attentive for further questions throughout the review process.\n\n\nReferences:\n\n[1] Wang, S., Yu, X., & Perdikaris, P. (2020). When and why PINNs fail to train: A neural tangent kernel perspective. arXiv preprint arXiv:2007.14527.\n\n[2] De Ryck, T., Bonnet, F., Mishra, S., & de Bézenac, E. (2024). An Operator Preconditioning Perspective on Training in Physics-Informed Machine Learning. ICLR 2024.", "claims_and_evidence": "The authors show that RTNN is constructed to satisfy the divergence-free form of systems of conservation laws. This is shown by straightforward calculation.\n\nIn numerical experiments, the authors illustrate that RTNN can be used in a collocated PDE-related loss minimization framework similar to Physics-Informed Neural Networks (PINNs) framework to solve various PDEs.", "methods_and_evaluation_criteria": "The evaluation criteria is relative L2 error against reference solution, which is reasonable.", "theoretical_claims": "The RTNN is constructed to theoretically satisfy the DSFT constraints, and the proof is straightforward. So numerically the fields satisfy the constraints to machine precision.", "experimental_designs_or_analyses": "The tested set of PDEs is significant that can be challenging for standard PINNs. However, it is not easy to see what kind of singular behaviors (sharp gradients) the solution has. Perhaps the authors can describe more in detail what kind of features are challenging to learn in these problems for existing methods like PINNs, and what contributes to the high errors for these methods (as opposed to RTNNs). There are some details in the current form of the manuscript, but it is hard to infer what the challenge is for the individual problems.", "relation_to_broader_scientific_literature": "Structure preservation is a major theme in classical numerical schemes, perhaps the authors can mention a few references (texts, survey papers).", "essential_references_not_discussed": "This type of embedding of conservation laws appear to be new, as far as I know.", "other_strengths_and_weaknesses": "The PDEs that appear in numerical examples are challenging PDEs for PINNs. The results show a clear advantage compared to standard PINNs.\n\nThe universal approximation theorem (Theorem 2.4) appears to be a very generic result, however, and does not shed a lot of light into why RTNNs have superior numerical performance.", "other_comments_or_suggestions": "I suggest the authors show a few of the solution plots (rather than the error plots only) so the readers can infer the difficulty in PINNs solution", "questions_for_authors": "Is there any numerical issues for approximating the solution that contains what is essentially a \"differentiation layer\"? Repeated differentiations can result in loss of significant digits, and can be a source of instability.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "r1vFlLh6kD", "reviewer": "Reviewer_G2wL", "summary": "In this paper, a method for learning divergence-free symmetric tensors is proposed.  As neural network architectures that maintain conservation laws, a typical approach is learning divergence-free vector fields by using the exterior calculus. The proposed method is another approach based on a different geometric structure. As far as I know, the use of this structure is certainly new. In numerical experiments, it is shown that the proposed method achieves higher accuracy than a method for learning divergence-free vector fields. I think that the proposed approach is very interesting and powerful.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for their time and feedback. Below we address some of the raised points:\n\n**Main Difference between NCL and RTNN:**  \nNCL enforces only the continuity (mass-conservation) equation, whereas our approach reformulates the conservation laws in flux form, constraining each flux tensor to be both divergence-free and symmetric. This flux representation allows us to extract primitive fields—density, momenta, and stress tensors—that are geometrically coupled and, by construction, satisfy both mass and momentum conservation. In fact, when the RTNN-derived fields are substituted into the flux-form momentum equation, the momentum conservation learning problem is reduced to an optimization problem on the stress tensors.\n\n**On NCL’s Performance in the Compressible Experiment:**  \nIn our compressible experiment, NCL underperformed compared to other approaches. In this case, we have variable density and an additional energy equation, a scenario that differs from the incompressible cases where NCL performed well (both in our experiments and in the original NCL paper).\n\nUpon investigation we observed that, in the compressible setting, the NCL architecture struggled to simultaneously satisfy the PDE residual and the initial condition, likely due to imbalances in the loss terms. To verify this hypothesis, we transformed the ansatz to the following form:\n\n$$\nF(t, x) = \\text{NCL}(t, x) - \\text{NCL}(0, x) + u_0(x)\n$$\n\nHere, \\\\( u_0(x) \\\\) represents the initial condition. This formulation enforces both the continuity and initial conditions by design, resulting in significantly improved accuracy and confirming our conjecture, and also supports the argument that hard constraint satisfaction leads to simpler optimization problems.\n\n**Related References:**  \nWe reviewed the reference by Liu *et al.* [1] and find it both interesting and relevant to our work, we thank the reviewer for bringing it to our attention. We will add it to the related work section of the camera-ready version.\n\nWe hope to have addressed the reviewer's question sufficiently and remain at their disposal throughout the review process.\n\n\n**References:**\n\n[1] Ning Liu, Yiming Fan, Xianyi Zeng, Milan Klöwer, Lu Zhang, Yue Yu, Harnessing the Power of Neural Operators with Automatically Encoded Conservation Laws, https://arxiv.org/abs/2312.11176", "claims_and_evidence": "The method is clearly explained with its theoretical background.", "methods_and_evaluation_criteria": "I believe that the proposed method is evaluated in an appropriate manner.", "theoretical_claims": "I checked the outline of the proofs, but I did not check them thoroughly.", "experimental_designs_or_analyses": "I believe that the experiments are designed appropriately; however, it seems strange that NCL performs worse than ordinary PINNs as seen in Table 1. In my understanding, similar to the proposed method, NCL is designed to satisfy conservation laws. However, the performance of one method (the proposed method) improves greatly while the other (NCL) declines, even though they are trying to maintain similar properties.", "relation_to_broader_scientific_literature": "This paper proposes a new structure for neural networks that preserves conservation laws. In previous papers, the property d d = 0 of the exterior derivative has been mainly used. The approach in this paper is different from those in the existing papers and has a potential to improve neural networks for physics.", "essential_references_not_discussed": "Perhaps discussing the following paper, which is an operator learning method that is based on a similar approach to CLN, is not essential but helpful for understanding the contributions of this paper.\n\nNing Liu, Yiming Fan, Xianyi Zeng, Milan Klöwer, Lu Zhang, Yue Yu, Harnessing the Power of Neural Operators with Automatically Encoded Conservation Laws, https://arxiv.org/abs/2312.11176", "other_strengths_and_weaknesses": "Overall, I think that this is a very interesting paper.", "other_comments_or_suggestions": "None.", "questions_for_authors": "I want to confirm the main difference between NCL and the proposed method. In my understanding, it is difficult for NCL to preserve multiple conservatoin laws simultaniously, but the proposed method can preserve them. Is this understanding correct?\n\nAlso, why did NCL perform so poorly in the numerical experiment?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "oNU1PRIqdE", "reviewer": "Reviewer_1n14", "summary": "This paper is within the setting of learning physical systems from PDE equations. The authors propose a neural network architecture that outputs a divergence-free symmetric tensor (DSFT). By adding this constraint, the neural network is guaranteed to enforce conservation of mass and momentum. In multiple experiments, the authors describe how to recover the solution from the DSFT output, after training the model on a loss function comprising PINNs loss, energy loss, and supervised L2 loss at some points. The authors compare their method to PINNs and Neural Conservation Laws (NCL), the latter of which is a similar method that only enforces conservation of mass. They also use their method to augment and compare to separable physics informed neural networks S-PINNs.\n\n# Post rebuttal response\nI raised my score from 2 to 3 after the authors' rebuttal (see my comment below).", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time reviewing the manuscript and their feedback that will contribute to improving the clarity of our work. Below we address the reviewer's concern:\n\n**Methods and Evaluation Criteria**:\n\n\nPINN and NCL were both trained using the full system (11)--(14), which includes the conservation of mass, momentum, and energy equations. For NCL , the mass conservation equation is satisfied by construction and cancels up to machine precision. \n\nFor RTNN, the continuity equation similarly cancels out and the energy equation (14) is also included in the loss. The only apparent difference is that RTNN replaces the explicit momentum residual with an isotropic loss term. However, this isotropic loss term naturally emerges when substituting the RTNN-predicted fields into the flux-form momentum equation. Imposing the zero-deviatoric condition on the stress tensor is numerically equivalent to enforcing the momentum equation directly in its expanded form (12--13). This equivalence can be evidenced following the reasoning below: \n\nConsider the compressible Euler system in (11)--(14):\n\nIn our RTNN parametrization, the network outputs a flux tensor \\\\(S\\\\) defined as:\n\n$$\nS = \\\\begin{pmatrix}\n\\\\rho & (\\\\rho u)^T \\\\\\\\\n\\\\rho u & \\\\rho u \\\\otimes u + \\\\sigma\n\\\\end{pmatrix},\n$$\n\nwith \\\\(\\\\rho = S_{0,0}\\\\) and \\\\(\\\\rho u = S_{0,1:n}\\\\) . The stress tensor is then defined by:\n\n$$\n\\\\sigma = S_{1:n,1:n} \\\\; - \\\\; \\\\rho\\\\, u \\\\otimes u.\n$$\n\nWe impose a zero-deviatoric constraint on \\\\(\\\\sigma\\\\), meaning its deviatoric part vanishes:\n\n$$\n\\\\sigma_{\\\\mathrm{dev}} \\\\;=\\\\; \\\\sigma \\\\; - \\\\; \\\\frac{1}{n}\\\\,\\\\operatorname{tr}(\\\\sigma)\\\\,I \\\\;\\\\triangleq  0.\n$$\n\nThis forces:\n\n$$\n\\\\sigma \\\\triangleq \\\\frac{1}{n}\\\\,\\\\operatorname{tr}(\\\\sigma)\\\\,I \\\\quad \\\\triangleq \\\\quad pI,\n$$\n\nwith the pressure defined as:\n\n$$\np =\\\\ \\\\frac{1}{n}\\\\,\\\\operatorname{tr}(\\\\sigma).\n$$\n\nSubstituting \\\\(\\\\sigma \\\\triangleq pI\\\\) back into the flux tensor, we obtain:\n\n$$\nS = \\\\begin{pmatrix}\n\\\\rho & (\\\\rho u)^T \\\\\\\\\n\\\\rho u & \\\\rho u \\\\otimes u + pI\n\\\\end{pmatrix}.\n$$\n\nOur architecture enforces that \\\\(S\\\\) is divergence-free:\n\n$$\n\\\\nabla \\\\cdot S = \\\\begin{pmatrix}\n\\\\partial_t \\\\rho + \\\\nabla \\\\cdot (\\\\rho u) \\\\\\\\[6pt]\n\\\\partial_t (\\\\rho u) + \\\\nabla \\\\cdot (\\\\rho u \\\\otimes u + pI)\n\\\\end{pmatrix} = 0.\n$$\n\nThe first row yields the continuity equation:\n\n$$\n\\\\partial_t \\\\rho + \\\\nabla \\\\cdot (\\\\rho u) = 0.\n$$\n\nThe second row yields:\n\n$$\n\\\\partial_t (\\\\rho u) + \\\\nabla \\\\cdot \\\\Bigl(\\\\rho u \\\\otimes u + pI\\\\Bigr) = 0,\n$$\n\nwhich when expanded yields the following equations:\n\n$$\n\\\\partial_t \\\\rho + \\\\nabla \\\\cdot (\\\\rho u, \\\\rho v) = 0 \\\\tag{11}\n$$\n\n$$\n\\\\partial_t (\\\\rho u) + \\\\nabla \\\\cdot (\\\\rho u^2, \\\\rho u\\\\,v) = -\\\\partial_x p \\\\tag{12}\n$$\n\n$$\n\\\\partial_t (\\\\rho v) + \\\\nabla \\\\cdot (\\\\rho u\\\\,v, \\\\rho v^2) = -\\\\partial_y p \\\\tag{13}\n$$\n\nThis is exactly the system in (11--13).\n\nA similar reasoning can be applied when adding viscous terms (or magnetic terms), for instance.\n\nFurthermore, to address the reviewer's concern numerically, we conducted an experiment in which the velocity and pressure fields extracted from a single RTNN were trained using the expanded momentum (12)--(13) equation rather than the simplified condition on the stress tensor. In this setting, the only difference among the experiments is the architecture. Although this process is more computationally expensive --- owing to the necessity of computing a fourth-order derivative --- the RTNN-trained fields still demonstrate similarly superior accuracy compared to those obtained with PINN and NCL. We summarize the obtained results below for the two unsteady Navier--Stokes problems: flow around a cylinder and the Beltrami flow for five different initialization seeds:\n\n| Method                                | Beltrami rL2 (Median ± IQR)         | Cylinder rL2 (Median ± IQR)         |\n|---------------------------------------|-------------------------------------|-------------------------------------|\n| RTNN (expanded momentum equation)     | 4.55e-04 ± 8.1e-05                  | 5.13e-03 ± 4.2e-04                  |\n| RTNN (stress tensor condition)        | 4.28e-04 ± 5.9e-05                  | 5.70e-03 ± 4.3e-04                  |\n| NCL                                   | 1.73e-03 ± 8.5e-05                  | 2.54e-02 ± 1.9e-03                  |\n| PINN                                  | 1.41e-03 ± 3.1e-04                  | 2.99e-02 ± 7.9e-04                  |\n\nWe acknowledge that the derivation of the condition on the stress tensors from the momentum equations could have been made clearer in the manuscript. We will appropriately make this link clear in the camera-ready version and provide the direct derivation in the appendix.\n\nWe hope we have addressed the reviewer's main concerns and remain at their disposal should further clarifications be needed.", "claims_and_evidence": "The method is presented in a straightforward way and makes sense that it would be an improvement over unconstrained PINNs and NCLs which only conserve mass. The experiments appear sensible and show a marked improvement over the selected baselines.", "methods_and_evaluation_criteria": "Yes, the evaluation criteria make sense. They compare to other common problems and show an order of magnitude improvement in each case. My expertise is not in PDEs, so I cannot attest to whether these problems are meaningful to that domain; this would be something I would defer to other reviewers on.", "theoretical_claims": "I did not check the proofs of the theoretical claims.", "experimental_designs_or_analyses": "I read through the design and setup of the experimental designs. I have a concern with how the competing methods are constructed. For example, in section 3.2, the loss function used to train the method is defined in terms of isotropic loss and energy conservation, rather than the PDE system described in (11)-(13). It seems that a different loss function is used for PINNs and NCL (line 247, column 2). Thus, it is unclear whether the performance increase is due to a different loss function or due to the novel architecture (the author's contribution).", "relation_to_broader_scientific_literature": "This method offers a novel architecture to provide hard constraints on conservation of mass and momentum. It builds on methods like NCL by hard-coding physical invariants into the architecture.", "essential_references_not_discussed": "I am not aware.", "other_strengths_and_weaknesses": "I think the idea to build constraints into the architecture makes sense and is a good approach. If my concern about the experimental section is addressed, I would consider raising my score.", "other_comments_or_suggestions": "N/a", "questions_for_authors": "Please see my comments under \"Methods and Evaluation Criteria\". While the results of your method appear strong, for the comparison to be fair the only difference between each comparator should be the architecture (i.e., PINNS=unconstrained, NCL=mass only, DSFT=mass and momentum). Please clarify whether this is the case in your experiments. If not, I request that you provide the results showing differences in architecture. Otherwise, it is unclear whether your method does better simply by using a specialized loss function for each problem (while NCL and PINNs do not).", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper proposes RTNNs, which can encodes the divergence-free constraints in neural networks within the PINN framework. The results show a clear advantage compared to standard PINNs. Some major concerns on experimental comparisons and evaluation metrics were successfully addressed by the authors' rebuttal.  These should be incorporated into the camera-ready version.
Decision: Accept (poster)