metrics, and examples, we expect to demonstrate the reliability and
versatility of the method.

###  4.1 Convergence analysis

We perform convergence tests for our method using a combination of analytical
results and numerical experiments. We break down the method into different
steps, including the Lie advection of the nematic field and the solution of
the surface Stokes equations.

Figure 4: Convergence test of different aspects of the method: A) Nematic sL
advection error with mesh refinement; B) Gaussian curvature error with mesh
refinement; C) Area dilation of the fluid velocity solution after each AL
pressure iteration; and D) Error of Stokes equations solution with mesh
refinement.

We start by testing the spatial convergence of the nematic semi-Lagrangian
advection on a unit sphere. To do this, we parameterize the unit sphere using
spherical coordinates, with polar angle $\theta\in[0,\pi]$ and azimuthal angle
$\varphi\in[0,2\pi)$, which results in the spherical vector bases,
$\boldsymbol{e}_{\theta}\coloneqq\partial_{\theta}/|\partial_{\theta}|$ and
$\boldsymbol{e}_{\varphi}\coloneqq\partial_{\varphi}/|\partial_{\varphi}|$. We
prescribe an analytical nematic field along the unit polar basis vector
$\boldsymbol{e}_{\theta}$, which creates two $+1$ defects at the poles of the
sphere. We also prescribe an advecting velocity $\boldsymbol{u}$, which is
defined using solenoidal vector spherical harmonics with $l=2$, $m=0$ (i.e.,
$\boldsymbol{u}=\boldsymbol{\phi}_{20}=(3/2)\sqrt{5/\pi}\sin\theta\cos\theta~{}\boldsymbol{e}_{\theta}$).
We then advect the nematic field for one time step using step size $\tau=0.1$
to measure the local error. We use the numerical result obtained from the
finest mesh tested as the ground truth, denoted as
$\overline{\boldsymbol{Q}}$, and compare results from coarser meshes using the
$L_{2}$ relative error norm,
$\varepsilon_{\boldsymbol{Q}}\coloneqq\|\boldsymbol{Q}-\overline{\boldsymbol{Q}}\|/\|\overline{\boldsymbol{Q}}\|$.
As shown in Fig. 4A, the error decreases approximately linearly with mesh
refinement, following $\varepsilon_{\boldsymbol{Q}}\propto h^{1.21}$, where
$h$ represents the average edge length of the mesh. Note that here we omitted
the unit-length normalization step in Alg. 1 and directly compared the
embedded $\mathbb{R}^{3\times 3}$ matrix representation of $\boldsymbol{Q}$.

The Stokes momentum equation involves the estimation of Gaussian curvature
$K$. To test this, we use a torus for which the analytical expression of
Gaussian curvature is known, with a value that varies from negative to
positive. We consider a torus with a radius of revolution $R$, and a radius of
cross-section $r$, for which the Gaussian curvature can be parametrized using
the polar coordinate as
$\overline{K}(\theta)={\cos\theta}/{(rR+r^{2}\cos\theta)}$. We calculate the
face-based Gaussian curvature $K$ using Eqn. (49) and compare it to the
analytical results using a relative $L_{2}$ error, denoted as
$\varepsilon_{K}\coloneqq\|K-\overline{K}\|/\|\overline{K}\|$, where
$\overline{K}(\theta)$ is computed at the center of the face. Fig. 4B
demonstrates that the error decreases linearly with mesh refinement, following
$\boldsymbol{\varepsilon}_{\boldsymbol{Q}}\propto h^{1.14}$.

We also evaluate the procedure of the AL that updates the pressure $p$ to
enforce the incompressibility constraint. This involves measuring the pressure
error norm that is defined by the the norm of area dilation,
$\varepsilon_{p}\coloneqq|\mathrm{div}~{}\boldsymbol{u}|$. Fig. 4C illustrates
that increasing the penalty coefficients $k$ results in faster convergence to
a divergence-free solution. However, this comes at the cost of making the
matrix inversion ill-conditioned, which increases the time per iteration.

Finally, we compare the numerical solution of the Stokes equations with its
analytical counterpart. We construct the forcing term $\boldsymbol{f}$ in the
momentum equation using the divergence-free spherical harmonics basis
function, which is an eigenfunction of the Stokes equations. This approach
allows us to derive an analytical solution $\bar{\boldsymbol{u}}$ for the
surface flow on a sphere. We chose $\boldsymbol{f}$ to be a mix of two
eigenmodes, $\boldsymbol{f}=\boldsymbol{\phi}_{21}+\boldsymbol{\phi}_{20}$,
and generated numerical solutions $\boldsymbol{u}$ across different mesh
resolutions. The error, as illustrated in Fig. 4D, decreases linearly with the
refinement of the mesh, following the relation
$\varepsilon_{\boldsymbol{u}}\coloneqq\|\boldsymbol{u}-\overline{\boldsymbol{u}}\|/\|\overline{\boldsymbol{u}}\|\propto
h^{1.34}$.

###  4.2 Computational time complexity

Figure 5: Computational time complexity of AL iterations for solving the
Stokes equations: the linear solve scales linearly with the system size $n$,
or the total number of faces in the mesh.

The Lie advection of the nematic field is explicit. The computation is also
local and parallelizable, which results in minimal computational cost. The
principal cost of the algorithm stems from the assembly of differential
matrices and the subsequent linear solve associated with the AL iteration. The
assembly of matrices, however, needs to be precomputed only once. Thus, the
predominant computation load is attributed to the Poisson solve A x = b for
the Stokes equations in each AL iteration. In this study, A is a symmetric
positive definite Laplacian matrix. Efficient solutions can be attained by
standard sparse linear algebra libraries. For smaller system sizes, we employ
the direct solve in Python’s Scipy.sparse library. As the system size
increases, an iterative method such as conjugate gradient method becomes more
pertinent. Empirically, as observed in Fig. 5, the computation time for the AL
iteration exhibits linear scaling with the system size $n$, where $n$ is the
total number of faces in the mesh.

These tests were conducted on a Dell 8940 Desktop equipped with an Intel CPU
i9-11900K operating at 3.50 GHz, 64 GB of RAM. In instances where a finer mesh
is employed, comprising approximately 10,000 vertices—akin to the examples
illustrated in the paper—the computation at each time step requires
approximately 2 minutes.

###  4.3 Nematodynamics on arbitrary surfaces

To validate our proposed methodology and demonstrate its versatility across a
spectrum of shapes with varying geometry and topology, we conducted a series
of numerical experiments. These experiments encompassed analytical shapes such
as spheres and genus-1 tori, as well as organic shapes characterized by
arbitrary curvature distributions. As previously mentioned in Secs. 2.2.1 and
3.6, the Stokes equations have Killing vector fields as their kernel. On
spheres and tori where such a Killing vector field exists, we carried out a
projection to remove the Killing component and focus on a specialized solution
with a minimal $L_{2}$ norm, $\|\boldsymbol{u}\|$. We present each example
under two distinct levels of activity, distinguished by disparate Péclet
numbers ($\mbox{\em Pe}=1$ and $\mbox{\em Pe}=10^{4}$) and effectively
showcasing both regular solutions and solutions in the regime of active
turbulence.

Figure 6:  Nematodynamics on a sphere under two levels of activity, as
captured by the Péclet number: A) low activity ($\mbox{\em Pe}=1$), and B)
high activity ($\mbox{\em Pe}=10^{4}$). The animated time series simulation is
available on Github: https://github.com/CunchengZhu/Riemannian-active-
nematics-2024/blob/main/sphere_compressed.mp4. From top to bottom, the time
plots show: the total numbers of $+1/2$ and $-1/2$ topological defects of the
nematic field, the $L_{2}$ norm squared $\|\boldsymbol{u}\|^{2}$ of the fluid
velocity, and the enstropy $\|\mathcal{K}\boldsymbol{u}\|^{2}$ of the surface
flow. Snapshots of the system are captured at the conclusion of each set of
simulations. The upper snapshots display the nematic field, color-mapped by
the local power density exerted by the nematic onto the fluid, expressed as
$P=-\langle\boldsymbol{u},\mathrm{div}^{\nabla}\boldsymbol{Q}\rangle$. The
lower snapshots visualize the velocity field, color-mapped by its magnitude
$|\boldsymbol{u}|$.

To characterize the nematodynamics, we focused on recording key metrics,
including the count of positive and negative defects (of charge $\pm 1/2$),
the $L_{2}$ norm squared of the flow velocity, $\|\boldsymbol{u}\|^{2}$, and
the entrosphy of the fluid, $\|\mathcal{K}\boldsymbol{u}\|^{2}$. The entrosphy
quantifies the magnitude of viscous dissipation and serves to assess the sharp
features and chaotic behavior of the flow field. For each case at each Péclet
number, we appended a visual snapshot of the nematodynamics simulation at the
end of each set of simulations. Taking Fig. 6 as an example, the top snapshot
displays the texture of the nematic field, with positive defects shown in red
and negative defects in blue. The bottom snapshot illustrates the streamlines
of the fluid flow field, color-mapped to represent velocity magnitude,
$|\boldsymbol{u}|$. The colormap on the nematic field represents the local
power density exerted by the active nematic on the fluid, expressed as
$P=-\langle\boldsymbol{u},\mathrm{div}^{\nabla}\boldsymbol{Q}\rangle$. Each
positive half defect is accompanied by a self-propelling fluid jet and a
vortex pair, and is surrounded by a region of positive power $P$. In contrast,
negative half defects exhibit 3-atic symmetry and generate an index-2 flow
field that is self-counteracting and remains static [68], contributing little
energy to the fluid on large scales.

In Fig. 6A, we examine dynamics on a sphere at Péclet number $\mbox{\em
Pe}=1$, corresponding to a low activity level. Here, we observe a fixed number
of four $+1/2$ defects exhibiting periodic orbiting motions–—a phenomenon
supported by previous empirical studies and simulations [11, 3]. In Fig. 6B,
we increase the activity level to $\mbox{\em Pe}=10^{4}$, leading to the
emergence of active turbulence. Active turbulence is characterized by unsteady
and chaotic dynamics driven by the spontaneous formation and annihilation of
topological defect pairs. Demonstrated by Fig. 6B and consistent with the
existing literature, the onset of turbulence is characterized by the
appearance of large bending deformation walls that separate aligned nematic
regions, the creation and annihilation of nematic defect pairs, and the
separation of positive defects from negative defects due to self-propulsion of
the positive [69]. As the total number of topological defects increases,
intricate flow patterns develop at smaller length scales. This observation is
supported by the fluctuating but steady enstrophy and a rapid decrease in the
overall velocity magnitude in the time series plots.

Figure 7: Nematodynamics on a torus at $\mbox{\em Pe}=1$ and $\mbox{\em
Pe}=10^{4}$. The animated time series simulations is available on Github:
https://github.com/CunchengZhu/Riemannian-active-
nematics-2024/blob/main/torus_compressed.mp4. See Fig. 6 for a detailed
caption.

As our method is able to accommodate different topologies, we expand our
exploration of nematodynamics to a genus one torus, characterized by a zero
Euler characteristic. The vanishing Euler characteristic allows for a defect-
free configuration at lower activity levels, as the nematic diffusion
gradually merges two pairs of positive and negative defects (cf. Fig. 7A). A
key observation is that the equilibrium defect-free nematic field on the torus
displays rotational symmetry around its vertical axis, while maintaining a
constant angle with the lines of toroidal coordinates. Consequently, the
resulting Stokes flow also exhibits this rotational symmetry, with the
toroidal midline acting as a shear layer. This flow pattern, in conjunction
with diffusion, stabilizes the nematic field, allowing the system to reach a
steady state. In Fig. 7B, when the Péclet number increases, the system
transitions into a state of active turbulence through the injection of energy
by topological defects, similar to what occurs on a sphere. Once again, we
note that enstropy reaches a plateau soon after the initial development of
defect bindings, while the velocity norm decreases as vortices form on
increasingly smaller length scales.

Figure 8: Nematodynamics on a topological sphere of varying curvature at
$\mbox{\em Pe}=1$ and $\mbox{\em Pe}=10^{4}$. The animated time series
simulations is available on Github: https://github.com/CunchengZhu/Riemannian-
active-nematics-2024/blob/main/spike_compressed.mp4. See Fig. 6 for a detailed
caption.

In addition to standard analytical shapes, we applied our method to a
topologically spherical membrane structure with an arbitrary curvature
distribution (cf. Fig. 8A). This geometry was generated through a simulation
based on the Helfrich bending energy under the influence of an osmotic shock
[70]. Even under low-activity conditions where nematic ordering relaxation
prevails, a relatively large total number of defects persists. Notably, these
defects tend to cluster in areas with like-sign curvature; that is, positive
defects gather at the protrusions with positive curvature, while negative
defects accumulate in negatively curved valleys. When two positive defects are
present at a protrusion, they create strong circulating flows at the tips.
These flows cause the defects to orbit periodically around the protrusions.
This interaction at the protrusion is isolated from the rest of the domain,
which remains almost motionless, thereby creating a limit cycle in the system.
The distinct behaviors of active nematics in areas of positive and negative
curvature suggests a possible detailed comparative study on active nematic
systems in hyperbolic and elliptic geometries as a followup work. More
broadly, the division of the flow field into areas of positive and negative
curvature has implications for understanding how mixing and coherent
structures are influenced by the shape of the domain. In the high activity
regime, there is a notable difference in flow velocities when comparing
conditions on a complex-shaped object to the sphere case. On a sphere, the
interaction of topological defects creates complex flow patterns on small
length scales, which tends to reduce the characteristic magnitude of the fluid
velocity. However, as seen in Fig. 8B, for the complex shape, the magnitude of
velocities remains relatively consistent whether in a steady state at low Pe
or in active turbulence at high Pe. As Pe increases, active turbulence
disrupts the stable, coherent structures that are present at low Pe, leading
to a state with more efficient mixing by the flow. Specifically, this results
in a decrease in velocity at the protrusions and an increase in velocity in
the valleys of the manifold.

##  5 Conclusion

In this study, we formulate a minimal model for the dynamics of an active
nematic fluid on a Riemannian 2-manifold. This model is cast based on a
coordinate-free differential-geometric language and on the theory of the
complex line bundle, which can naturally be generalized to curved spaces. We
use the complex line to establish the 2D $k$-atic equivalence and model
nematic advection through generalized Lie derivatives, laying the groundwork
for extending the model to $k$-atic hydrodynamics. The Levi-Civita connection
and its curvature are introduced within the framework of the complex line
bundle. We show the coupling of the fluid’s viscous stress with the curvature
through the covariant Dolbeault $\bar{\partial}^{\nabla}$ operator.
Numerically, our method evolves the active nematic system on a triangular mesh
with arbitrary shape and topology. This is done by adhering to the discrete
analogue of the continuous theory. We construct the Hermitian Bochner
Laplacian matrix based on the variational principle and develop a generalized
semi-Lagrangian scheme according to the Lie advection. This method maintains
robustness and efficiency across various flow regimes, from low to high
activity levels.

In future work, we anticipate detailed studies that examine the relationship
between the complex number representation and the well-established
$\boldsymbol{Q}$-tensor theory, including their generalizations to $k$-atic
fields. We also acknowledge that our minimal Riemannian nematodynamics model
does not account for the possible influence of extrinsic geometry on nematic
and fluid dynamics, which may be significant in various applications. In
subsequent work that will be detailed elsewhere, we explore the surface Stokes
flow of evolving manifolds, where fluid flow induces extrinsic deformations of
the surface. Depending on the specific application, it may also be interesting
to include factors such as membrane bending rigidity, interactions with bulk
fluid, and reorientation of nematic molecules along the extrinsic principal
curvatures.

In summary, our integration of geometric language and complex manifold theory
into the study of active nematics sets the stage for further theoretical and
computational developments in this area. With further post-processing and
analysis, we expect that the computational tool developed here will also be
instrumental for specialists seeking deeper understanding of active nematics
influenced by geometric aspects in specific systems. We hope that the
incorporation of realistic geometric and topological features in the model can
further narrow the gap between theoretical predictions and experimental
findings.

## Acknowledgements

C.Z. thanks M. Firouznia and Y. Chen for their helpful references and many
discussions on various topics covered in this study. D.S. acknowledges funding
from National Science Foundation Grant No. DMS-2153520. A.C. acknowledges
funding from National Science Foundation CAREER Award 2239062. Additional
support was provided by SideFX software.

##  Appendix A Derivation of $\bar{\partial}^{\nabla*}\bar{\partial}^{\nabla}$
using complex differential form

In this section, we derive the viscosity Laplacian,
$-2\bar{\partial}^{\nabla*}\bar{\partial}^{\nabla}=\Delta+K:\Gamma(TM)\rightarrow\Gamma(TM)$,
through the use of complex differential forms, assuming a basic knowledge of
exterior calculus. Consider a Hermitian complex line bundle $L$ over $M$ with
connection $\nabla\colon\Gamma(L)\to\Omega^{1}(M;L)$, where $\Omega^{k}(M;L)$
represents the space of $L$-valued $k$-forms. The curvature 2-form for the
connection, $R^{\nabla}\in\Omega^{2}(M;\mathbb{C})$, is defined by

$\displaystyle
R^{\nabla}\llbracket\boldsymbol{v},\boldsymbol{w}\rrbracket\boldsymbol{\psi}\coloneqq(d^{\nabla}\nabla\boldsymbol{\psi})\llbracket\boldsymbol{v},\boldsymbol{w}\rrbracket=\nabla_{\boldsymbol{v}}\nabla_{\boldsymbol{w}}\boldsymbol{\psi}-\nabla_{\boldsymbol{w}}\nabla_{\boldsymbol{v}}\boldsymbol{\psi}-\nabla_{[\boldsymbol{v},\boldsymbol{w}]}\boldsymbol{\psi}.$
(A.1)

Note that when $L=TM$ and $\nabla$ is the Levi-Civita connection, $R^{\nabla}$
correlates directly with the Gaussian curvature $K$ as
$R^{\nabla}=-\boldsymbol{\mathrm{i}}K\star 1$.

Note that for a 1-form $\alpha\in\Omega^{1}(M;L)$ on a 2-manifold, the
relation
$\alpha\llbracket\boldsymbol{\mathrm{i}}\boldsymbol{v}\rrbracket=-\star\alpha\llbracket
v\rrbracket$ holds, where $\star$ denotes the Hodge star operator. The anti-
holomorphic derivative, $\bar{\partial}^{\nabla}\boldsymbol{\psi}$,
$\boldsymbol{\psi}\in\Gamma(L)$, induced by the connection $\nabla$, can be
expressed as

$\displaystyle\bar{\partial}^{\nabla}\boldsymbol{\psi}\coloneqq\frac{1}{2}(\nabla\boldsymbol{\psi}+\boldsymbol{\mathrm{i}}\nabla_{\boldsymbol{\mathrm{i}}(\cdot)}\boldsymbol{\psi})=\frac{1}{2}(\nabla\boldsymbol{\psi}-\boldsymbol{\mathrm{i}}\star\nabla\boldsymbol{\psi}).$
(A.2)

The anti-holomoprhic derivative allows us to construct the anti-holomorphic
energy as follows:

$\displaystyle\begin{split}2\langle\bar{\partial}^{\nabla}\boldsymbol{\psi},\bar{\partial}^{\nabla}\boldsymbol{\psi}\rangle\star
1&=\frac{1}{2}\langle\nabla\boldsymbol{\psi}-\boldsymbol{\mathrm{i}}\star\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}-\boldsymbol{\mathrm{i}}\star\nabla\boldsymbol{\psi}\rangle\star
1\\\
&=\frac{1}{2}\left(\langle\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}\rangle\star
1+\langle\boldsymbol{\mathrm{i}}\star\nabla\boldsymbol{\psi},\boldsymbol{\mathrm{i}}\star\nabla\boldsymbol{\psi}\rangle\star
1-2\langle\boldsymbol{\mathrm{i}}\star\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}\rangle\star
1\right)\\\
&=\langle\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}\rangle\star
1+\langle\nabla\boldsymbol{\psi}{\wedge}\boldsymbol{\mathrm{i}}\nabla\boldsymbol{\psi}\rangle\\\
&=\langle\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}\rangle\star
1+d\langle\boldsymbol{\psi}\wedge\boldsymbol{\mathrm{i}}\nabla\boldsymbol{\psi}\rangle-\langle\boldsymbol{\psi}\wedge\boldsymbol{\mathrm{i}}d^{\nabla}\nabla\boldsymbol{\psi}\rangle\\\
&=\langle\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}\rangle\star
1-\langle\boldsymbol{\psi},\star\boldsymbol{\mathrm{i}}d^{\nabla}\nabla\boldsymbol{\psi}\rangle=\langle\nabla\boldsymbol{\psi},\nabla\boldsymbol{\psi}\rangle\star
1-\langle\boldsymbol{\psi},\boldsymbol{\mathrm{i}}\star
R^{\nabla}\boldsymbol{\psi}\rangle.\end{split}$ (A.3)

By taking a variation with respect to $\boldsymbol{\psi}$ on the final
expression, we obtain the operator
$-2\bar{\partial}^{\nabla*}\bar{\partial}^{\nabla}=\Delta+\boldsymbol{\mathrm{i}}\star
R^{\nabla}:\Gamma(L)\rightarrow\Gamma(L)$. When the complex line bundle is
specialized to the tangent bundle, $L=TM$, the operator is the viscosity
Laplacian,
$2\bar{\partial}^{\nabla*}\bar{\partial}^{\nabla}=\Delta+K:\Gamma(TM)\rightarrow\Gamma(TM)$.

## References

  * Alert et al. [2022] R. Alert, J. Casademunt, J.-F. Joanny,  Active turbulence,  Annual Review of Condensed Matter Physics 13 (2022) 143–170.
  * Doostmohammadi et al. [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, F. Sagués,  Active nematics,  Nature Communications 9 (2018) 3246.
  * Firouznia and Saintillan [2024] M. Firouznia, D. Saintillan,  Self-organized dynamics of a viscous drop with interfacial nematic activity,  arXiv preprint arXiv:2404.11729 (2024).
  * Apaza and Sandoval [2018] L. Apaza, M. Sandoval,  Active matter on Riemannian manifolds,  Soft Matter 14 (2018) 9928–9936.
  * Mietke et al. [2019] A. Mietke, F. Jülicher, I. F. Sbalzarini,  Self-organized shape dynamics of active surfaces,  Proceedings of the National Academy of Sciences 116 (2019) 29–34.
  * Hoffmann et al. [2022] L. A. Hoffmann, L. N. Carenza, J. Eckert, L. Giomi,  Theory of defect-mediated morphogenesis,  Science Advances 8 (2022) eabk2712.
  * Bell et al. [2022] S. Bell, S.-Z. Lin, J.-F. Rupprecht, J. Prost,  Active nematic flows over curved surfaces,  Physical Review Letters 129 (2022) 118001.
  * Al-Izzi and Morris [2021] S. C. Al-Izzi, R. G. Morris,  Active flows and deformable surfaces in development,  Seminars in Cell & Developmental Biology 120 (2021) 44–52.
  * Al-Izzi and Morris [2023] S. C. Al-Izzi, R. G. Morris,  Morphodynamics of active nematic fluid surfaces,  Journal of Fluid Mechanics 957 (2023) A4.
  * Zhang et al. [2016] R. Zhang, Y. Zhou, M. Rahimi, J. J. De Pablo,  Dynamic structure of active nematic shells,  Nature Communications 7 (2016) 13483.
  * Keber et al. [2014] F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, A. R. Bausch,  Topology and dynamics of active nematic vesicles,  Science 345 (2014) 1135–1139.
  * Napoli and Vergori [2012] G. Napoli, L. Vergori,  Extrinsic curvature effects on nematic shells,  Physical Review Letters 108 (2012) 207803.
  * Napoli and Vergori [2016] G. Napoli, L. Vergori,  Hydrodynamic theory for nematic shells: The interplay among curvature, flow, and alignment,  Physical Review E 94 (2016) 020701.
  * Sonnet and Virga [2012] A. M. Sonnet, E. G. Virga, Dissipative Ordered Fluids: Theories for Liquid Crystals, Springer Science & Business Media, 2012\.
  * Salbreux et al. [2009] G. Salbreux, J. Prost, J.-F. Joanny,  Hydrodynamics of cellular cortical flows and the formation of contractile rings,  Physical Review Letters 103 (2009) 058102.
  * Giomi et al. [2022] L. Giomi, J. Toner, N. Sarkar,  Hydrodynamic theory of p-atic liquid crystals,  Physical Review E 106 (2022) 024701.
  * Doi and Edwards [1988] M. Doi, S. F. Edwards, The Theory of Polymer Dynamics, Clarendon Press, 1988\.
  * Bird et al. [1988] R. B. Bird, O. Hassager, R. C. Armstrong, C. F. Curtiss, Dynamics of Polymeric Liquids: Kinetic Theory, Wiley Interscience, 1988.
  * Beris and Edwards [1994] A. N. Beris, B. J. Edwards, Thermodynamics of Flowing Systems With Internal Microstructure, Oxford University Press, 1994\.
  * Saintillan and Shelley [2013] D. Saintillan, M. J. Shelley,  Active suspensions and their nonlinear models,  Comptes Rendus Physique 14 (2013) 497–517.
  * Morozov and Spagnolie [2015] A. Morozov, S. E. Spagnolie,  Introduction to complex fluids,  in: S. E. Spagnolie (Ed.), Complex Fluids in Biological Systems: Experiment, Theory, and Computation, Springer, 2015, pp. 3–52.
  * Gilbert and Vanneste [2023] A. D. Gilbert, J. Vanneste,  A geometric look at momentum flux and stress in fluid mechanics,  Journal of Nonlinear Science 33 (2023) 31.
  * Wensch and Naumann [2012] J. Wensch, A. Naumann,  Semi-Lagrangian discretization of the upper-convective derivative in Non-Newtonian fluid flow,  Proceedings in Applied Mathematics and Mechanics 12 (2012) 765–766.
  * Nabizadeh et al. [2022] M. S. Nabizadeh, S. Wang, R. Ramamoorthi, A. Chern,  Covector fluids,  ACM Transactions on Graphics 41 (2022) 1–16.
  * Mullen et al. [2011] P. Mullen, A. McKenzie, D. Pavlov, L. Durant, Y. Tong, E. Kanso, J. E. Marsden, M. Desbrun,  Discrete lie advection of differential forms,  Foundations of Computational Mathematics 11 (2011) 131–149.
  * Pismen [1999] L. M. Pismen, Vortices in Nonlinear Fields: From Liquid Crystals to Superfluids, From Non-Equilibrium Patterns to Cosmic Strings, Oxford University Press, 1999.
  * Chan et al. [2017] C. H. Chan, M. Czubak, M. M. Disconzi,  The formulation of the Navier–Stokes equations on Riemannian manifolds,  Journal of Geometry and Physics 121 (2017) 335.
  * Samavaki and Tuomela [2020] M. Samavaki, J. Tuomela,  Navier–Stokes equations on Riemannian manifolds,  Journal of Geometry and Physics 148 (2020) 103543.
  * Rank and Voigt [2021] M. Rank, A. Voigt,  Active flows on curved surfaces,  Physics of Fluids 33 (2021) 072110.
  * Brandner et al. [2022] P. Brandner, T. Jankuhn, S. Praetorius, A. Reusken, A. Voigt,  Finite element discretization methods for velocity-pressure and stream function formulations of surface Stokes equations,  SIAM Journal on Scientific Computing 44 (2022) A1807.
  * Reuther et al. [2020] S. Reuther, I. Nitschke, A. Voigt,  A numerical approach for fluid deformable surfaces,  Journal of Fluid Mechanics 900 (2020) R8.
  * Zhang et al. [2016] R. Zhang, T. Roberts, I. S. Aranson, J. J. De Pablo,  Lattice Boltzmann simulation of asymmetric flow in nematic liquid crystals with finite anchoring,  Journal of Chemical Physics 144 (2016) 084905.
  * Nitschke et al. [2021] I. Nitschke, S. Reuther, A. Voigt,  Vorticity‐stream function approaches are inappropriate to solve the surface Navier‐Stokes equation on a torus,  Proceedings in Applied Mathematics and Mechanics 20 (2021) e202000006.
  * Yin et al. [2023] H. Yin, M. S. Nabizadeh, B. Wu, S. Wang, A. Chern,  Fluid cohomology,  ACM Transactions on Graphics 42 (2023) 126.
  * Reusken [2020] A. Reusken,  Stream function formulation of surface Stokes equations,  IMA Journal of Numerical Analysis 40 (2020) 109.
  * Pearce et al. [2019] D. Pearce, P. W. Ellis, A. Fernandez-Nieves, L. Giomi,  Geometrical control of active turbulence in curved topographies,  Physical Review Letters 122 (2019) 168002.
  * Torres-Sánchez et al. [2019] A. Torres-Sánchez, D. Millán, M. Arroyo,  Modelling fluid deformable surfaces with an emphasis on biological interfaces,  Journal of Fluid Mechanics 872 (2019) 218–271.
  * Gross et al. [2018] S. Gross, T. Jankuhn, M. A. Olshanskii, A. Reusken,  A trace finite element method for vector-Laplacians on surfaces,  SIAM journal on numerical analysis 56 (2018) 2406–2429.
  * Ben-Chen et al. [2010] M. Ben-Chen, A. Butscher, J. Solomon, L. Guibas,  On discrete Killing vector fields and patterns on surfaces,  Computer Graphics Forum 29 (2010) 1701–1711.
  * de Goes et al. [2014] F. de Goes, B. Liu, M. Budninskiy, Y. Tong, M. Desbrun,  Discrete 2-tensor fields on triangulations,  Computer Graphics Forum 33 (2014) 13–24.
  * De Goes et al. [2020] F. De Goes, A. Butts, M. Desbrun,  Discrete differential operators on polygonal meshes,  ACM Transactions on Graphics (TOG) 39 (2020) 110–1.
  * Knöppel et al. [2013] F. Knöppel, K. Crane, U. Pinkall, P. Schröder,  Globally optimal direction fields,  ACM Transactions on Graphics 32 (2013) 1–10.
  * Sageman-Furnas et al. [2019] A. O. Sageman-Furnas, A. Chern, M. Ben-Chen, A. Vaxman,  Chebyshev nets from commuting PolyVector fields,  ACM Transactions on Graphics 38 (2019) 1–16.
  * De Goes et al. [2016] F. De Goes, M. Desbrun, Y. Tong,  Vector field processing on triangle meshes,  in: ACM SIGGRAPH 2016 Courses, 2016, pp. 1–49.
  * Schutz [1980] B. F. Schutz, Geometrical Methods of Mathematical Physics, Cambridge university press, 1980\.
  * Kanso et al. [2007] E. Kanso, M. Arroyo, Y. Tong, A. Yavari, J. G. Marsden, M. Desbrun,  On the geometric character of stress in continuum mechanics,  Zeitschrift für angewandte Mathematik und Physik 58 (2007) 843–856.
  * Needham [2021] T. Needham, Visual Differential Geometry and Forms: A Mathematical Drama in Five Acts, Princeton University Press, 2021.
  * De Gennes and Prost [1993] P.-G. De Gennes, J. Prost, The Physics of Liquid Crystals, 83, Oxford University Press, 1993.
  * Frank [1958] F. C. Frank,  I. Liquid crystals. On the theory of liquid crystals,  Discussions of the Faraday Society 25 (1958) 19–28.
  * Needham [2023] T. Needham, Visual Complex Analysis, Oxford University Press, 2023.
  * Oldroyd [1950] J. G. Oldroyd,  On the formulation of rheological equations of state,  Proceedings of the Royal Society of London A 200 (1950) 523–541.
  * Doi [2013] M. Doi, Soft Matter Physics, Oxford University Press, 2013.
  * Doi [2011] M. Doi,  Onsager’s variational principle in soft matter,  Journal of Physics: Condensed Matter 23 (2011) 284118.
  * Arroyo et al. [2018] M. Arroyo, N. Walani, A. Torres-Sánchez, D. Kaurin,  Onsager’s variational principle in soft matter: Introduction and application to the dynamics of adsorption of proteins onto fluid membranes,  in: D. J. Steigmann (Ed.), The Role of Mechanics in the Study of Lipid Bilayers, Springer International Publishing, 2018, pp. 287–332.
  * Mirza et al. [2023] W. Mirza, A. Torres-Sánchez, G. Vilanova, M. Arroyo,  Variational formulation of active nematics: theory and simulation,  arXiv preprint arXiv:2306.01515 (2023).
  * Arroyo and DeSimone [2009] M. Arroyo, A. DeSimone,  Relaxation dynamics of fluid membranes,  Physical Review E 79 (2009) 031915.
  * Saintillan [2018] D. Saintillan,  Rheology of active fluids,  Annual Review of Fluid Mechanics 50 (2018) 563–592.
  * Gao et al. [2017] T. Gao, M. D. Betterton, A.-S. Jhang, M. J. Shelley,  Analytical structure, dynamics, and coarse graining of a kinetic model of an active fluid,  Physical Review Fluids 2 (2017) 093302.
  * Ericksen [1991] J. L. Ericksen,  Liquid crystals with variable degree of orientation,  Archive for Rational Mechanics and Analysis 113 (1991) 97.
  * Virga [2018] E. G. Virga, Variational Theories for Liquid Crystals, CRC Press, 2018.
  * Leslie [1968] F. M. Leslie,  Some constitutive equations for liquid crystals,  Archive for Rational Mechanics and Analysis 28 (1968) 265–283.
  * Canevari and Segatti [2018] G. Canevari, A. Segatti,  Defects in nematic shells: a $\gamma$-convergence discrete-to-continuum approach,  Archive for Rational Mechanics and Analysis 229 (2018) 125–186.
  * Bretherton [1962] F. P. Bretherton,  The motion of rigid particles in a shear flow at low Reynolds number,  Journal of Fluid Mechanics 14 (1962) 284–304.
  * Thijssen [2020] K. Thijssen,  Active nematics with anisotropic friction: the decisive role of the flow aligning parameter,  Soft Matter 16 (2020) 2065.
  * Jeffery [1922] G. B. Jeffery,  The motion of ellipsoidal particles immersed in a viscous fluid,  Proceedings of the Royal Society of London A 102 (1922) 161–179.
  * Hundsdorfer and Verwer [2003] W. Hundsdorfer, J. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Berlin Heidelberg, 2003.
  * Quarteroni and Valli [2008] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer Science & Business Media, 2008\.
  * Giomi et al. [2014] L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, M. Cristina Marchetti,  Defect dynamics in active nematics,  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372 (2014) 20130365.
  * Thampi and Yeomans [2016] S. Thampi, J. Yeomans,  Active turbulence in active nematics,  The European Physical Journal Special Topics 225 (2016) 651–662.
  * Zhu et al. [2022] C. Zhu, C. Lee, P. Rangamani,  Mem3DG: Modeling membrane mechanochemical dynamics in 3D using discrete differential geometry.,  Biophysical Reports 2 (2022) 100062–100062.


# Neural Networks with Inputs Based on Domain of Dependence and A Converging
Sequence for Solving Conservation Laws, Part I: 1D Riemann Problems

Haoxiang Huang11footnotemark: 1, Yingjie Liu22footnotemark: 2 and Vigor
Yang33footnotemark: 3

###### Abstract

Recent research works for solving partial differential equations (PDEs) with
deep neural networks (DNNs) have demonstrated that spatiotemporal function
approximators defined by auto-differentiation are effective for approximating
nonlinear problems, e.g. the Burger’s equation, heat conduction equations,
Allen-Cahn and other reaction-diffusion equations, and Navier-Stokes equation.
Meanwhile, researchers apply automatic differentiation in physics-informed
neural network (PINN) to solve nonlinear hyperbolic systems based on
conservation laws with highly discontinuous transition, such as Riemann
problem, by inverse problem formulation in data-driven approach. However, it
remains a challenge for forward methods using DNNs without knowing part of the
solution to resolve discontinuities in nonlinear conservation laws. In this
study, we incorporate 1st order numerical schemes into DNNs to set up the loss
functional approximator instead of auto-differentiation from traditional deep
learning framework, e.g. TensorFlow package, which improves the effectiveness
of capturing discontinuities in Riemann problems. In particular, the 2-Coarse-
Grid neural network (2CGNN) and 2-Diffusion-Coefficient neural network (2DCNN)
are introduced in this work. We use $2$ solutions of a conservation law from a
converging sequence, computed from a low-cost numerical scheme, and in a
domain of dependence of a space-time grid point as the input for a neural
network to predict its high-fidelity solution at the grid point. Despite
smeared input solutions, they output sharp approximations to solutions
containing shocks and contacts and are efficient to use once trained.

††footnotetext: Key words: 2-Coarse-Grid neural network; 2-Diffusion-
Coefficient neural network; physics informed machine learning; Riemann
problem; multi-fidelity optimization.††footnotetext: 1\. E-mail:
<EMAIL_ADDRESS>Woodruff School of Mechanical Engineering, Georgia
Institute of Technology, Atlanta, GA 30332, USA.††footnotetext: 2\. E-mail:
<EMAIL_ADDRESS>School of Mathematics, Georgia Institute of
Technology, Atlanta, GA 30332, USA. Research is supported in part by NSF grant
DMS-1522585.††footnotetext: 3\. Email<EMAIL_ADDRESS>Daniel
Guggenheim School of Aerospace Engineering, Georgia Institute of Technology,
Atlanta, GA 30332, USA.

##  1 Introduction

There has been a lot of research works on numerical methods for solving
conservation laws whose solutions may contain shocks and contact
discontinuities, such as the Godunov scheme [1], MUSCL scheme [2, 3], ENO [4,
5] and WENO [6, 7] schemes, hierarchical reconstruction [8, 9, 10], and many
others. Numerical techniques are important for studying high-speed aerodynamic
flows which plays a substantial role in aircraft designs, combustion problems
and astronomy physics [11, 12, 13, 14]. However, the development of machine
learning techniques for solving hyperbolic conservation laws are still in
early stage.

In the past decade, data driven modeling in machine learning has been widely
developed for multiple scientific disciplines, including image processing,
biomedical applications and engineering design optimization [15, 16, 17, 18,
19, 20, 21]. Especially, advances in computational resources, e.g. graphics
processing unit (GPU) and tensor processing unit (TPU), accelerate training
speed in deep learning frameworks, e.g. TensorFlow, PyTorch, for computer
vision, natural language processing and other important scientific disciplines
[22, 23, 24].

Recently, Raissi et al. [25] employed deep neural networks with physics
informed conditions, e.g. initial conditions, boundary conditions and
functional forms of partial differential equations (PDEs) by applying chain
rules for differentiating compositions of functions utilizing automatic
differentiation [26]. This method has achieved many successes in data driven
methods for identifying nonlinear PDEs or predicting turbulent mixing, vortex
induced vibration (VIV) with given governing PDEs, e.g., the Navier-Stokes
equations [27, 28, 29, 25, 30, 31]. Meanwhile, PINN has demonstrated that it
can solve more physical problems from reality, e.g., Lattice Boltzmann
Equation with the Bhatnagar-Gross-Krook collision, hypersonic flow, heat
transer, nano-optics, electro-convection and so on [32, 33, 34, 35, 36, 37].
Besides PINN developed based on DNN, recent researchers have discovered that
physics-informed convolutional neural network (CNN) are capable of resolving
physical problems with irregular boundary values [38]. Meanwhile, to explore
the possibility of learning temporal stages in the future while data may not
be available, Psaros et al. [39] extended PINN to the meta-learning framework.
In Nguyen and Tsai et al. [40], a neural network is developed for solving
second order wave equations by using the solution and its gradient computed
from a low order scheme as input to predict the higher order solution.

To examine the possibility of applying PINN to solving the Euler equations,
Mao et al. [41] designed forward and inverse methods by selecting clustered
training data points in the discontinuity region along with randomly selected
data points in other spatio-temporal positions, which enhance the accuracy of
predictions for the Lax and Sod problems. Another study by Michoski et al.
[42] compares solutions computed by finite volume methods and results from
predictions of PINN. In their PINN algorithm design, diffusion terms are
introduced in functional forms of Euler equations inside the neural network to
smooth out discontinuities, though there is no detailed discussion of the
functional forms. We are still wondering whether automatic differentiation
along works in defining functional forms of Euler equations inside the neural
network for solving the equations, without knowing any part of the solution.
Is there any efficient methodology using neural networks to compute the Euler
equations with results comparable to those from high resolution non-
oscillatory schemes? We would like to refer to prior works on locating shocks
or contacts from a numerical solution or smeared solution [43, 44].

In our study, we first apply automatic differentiation to approximate Euler
equations (with or without artificial diffusion) in the cost functional form
of PINN which doesn’t seem to produce good approximation of the solutions. To
resolve this issue, 1st order numerical schemes are used to approximate the
Euler equations in the cost functional form of PINN. They produce smeared
approximations to the solutions comparable to those obtained from
corresponding schemes, and are referred to as numerical PINN (N-PINN)
throughout the paper. To improve the resolution of N-PINN, scattered data
points are randomly selected as the additional training data in the vicinities
of discontinuities of the solution in the training process. This greatly
increases the resolution of N-PINN. Finally, we introduce our best neural
network methods which are efficient once trained and are able to produce sharp
results competitive to those from high resolution numerical schemes. In order
for the neural network to predict the solution at a space-time grid point, we
first compute two coarse grid 1st order numerical solutions and then use them
in a carefully selected small domain of dependence of the grid point as input.
Alternatively, instead of using two grids, we can use a low-cost scheme on one
grid to compute a conservation law perturbed with two diffusion coefficients,
and use the solutions in a carefully selected small domain of dependence of
the grid point as input. We find that the 2-coarse-grid neural network and
2-Diffusion-Coefficient neural network can reduce the error of prediction much
better than that from N-PINN even with selected partial data and at the same
time reduce training time substantially, which may save time for engineering
design optimization.

##  2 Problem Setup

Multi-layer neural networks [45] has found their application in solving PDEs,
e.g., through PINN which uses automatic differentiation (auto-differentiation)
to set up functional forms of PDEs [25]. However, the Euler equations remain a
challenge for this approach because of the complexity of the system and
discontinuities occurring in the solution. We first set up PDE functional
forms inside the neural network using 1st order numerical schemes, including
the Rusanov scheme, Leapfrog scheme with diffusion terms, and Leapfrog scheme
with nonlinearly weighted diffusion terms to approximate the solution.

In this study, we start with 1-D conservation laws for our proposed methods. A
scalar conservation law is defined as

$\frac{\partial U}{\partial t}+\frac{\partial f(U)}{\partial
x}=0,x\in\Omega\subset\mathcal{R},t\in[0,T],$ (1)

and the 1-D Euler equations for an ideal gas are

$\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\\ \rho v\\\
E\end{array}\right)+\frac{\partial}{\partial x}\left(\begin{array}[]{c}\rho
v\\\ \rho v^{2}\\\
v(E+p)\end{array}\right)=0,\;x\in\Omega\subset\mathcal{R},\;t\in[0,T]~{},$ (2)

where $\rho$, $u$ and $p$ are density, velocity and pressure respectively,
$\Omega$ is an interval,

$E=\frac{p}{\gamma-1}+\frac{1}{2}\rho v^{2}~{},$ (3)

$\gamma$ is the specific heat ratio and its value is $1.4$ throughout all the
cases studied in the paper.

To approximate the solution of the Euler equations by DNN, we use physics
informed dataset, e.g. initial conditions, boundary conditions and training
data from $\rho$, $v$, and $p$ to allow DNN to infer all states of interest in
the spatiotemporal domain $(x,t)$. We design a DNN with a spatiotemporal point
as its input and the prime variables ($\rho$, $v$, and $p$) of the predicted
solution at the point as its output. In order to define the loss function of
the DNN, we first partition the space-time domain with a rectangular grid,
then consider a $1^{st}$ order numerical scheme expressed in terms of the
prime variables at grid points (through the output of DNN) approximating each
equation of the Euler equations. The loss function is the summation of two
parts. The first part is the loss of the physics-informed dataset from initial
conditions, boundary conditions and training data (if available), while the
second part is the sum of square of the residue error from the 1st order
scheme at each grid point for each equation of the Euler system. The weights
and biases of the DNN can be learned from minimizing the loss function in this
PINN. The loss function is defined as follows.

$Err_{total}=Err_{1}+Err_{2},$ (4)

where $Err_{1}$ is the loss of physics-informed data defined as

$Err_{1}=\frac{1}{N_{1}}\sum^{N_{1}}_{i=1}[u(t^{i},x^{i})-u^{i}]^{2},$ (5)

superscript $i$ indicates a spatiotemporal point where the value of a prime
variable $u$ is available as $u^{i}$, $u(t^{i},x^{i})$ is the corresponding
DNN-predicted value at the point, $Err_{2}$ measures the residue error
(denoted as the $1^{st}$ order scheme functional form) defined as

$Err_{2}=\frac{1}{N_{2}}\sum^{N_{2}}_{j=1}|L(t^{j},x^{j})|^{2},$ (6)

superscript $j$ indicates a spatiotemporal grid point, $L(t^{j},x^{j})$ is the
residue of the scheme at the grid point (in which the solution values inside
the scheme are substituted by corresponding DNN-predicted values) and
$|\cdot|$ is the Euclidean norm. We denote
$\\{t^{i},x^{i},u^{i}\\}_{i=1}^{N_{1}}$ as the physics-informed dataset which
contains initial and boundary data (known as the forward problem) and possibly
additional training data (known as the inverse problem.)

##  3 Numerical Physics-Informed Neural Network (N-PINN)

The original PINN defines the functional form by auto-differentiation.
However, auto-differentiation seems to have difficulty in capturing
discontinuous solutions for a highly nonlinear system such as the Euler
system. We use PINN to predict the solution of the Euler system without any
training data point in the interior of the spatiotemporal domain. The
schematic below depicts the training procedure for PINN.

Figure 1: Schematic of physics informed neural network (PINN) with PDE
functional form (residue) approximated by auto-differentiation.

In this section, We define another PDE functional form for the Euler equations
by using a $1^{st}$ order numerical scheme, and the schematic describing the
training procedure is modified as followed:

Figure 2: Schematic of N-PINN with PDE functional form (residue) approximated
by a $1^{st}$ order numerical scheme.

The numerical scheme used in the functional form of N-PINN is either the
Rusanov scheme [46] or a Leapfrog-type scheme (applied to the Euler equations
with arificial diffusion terms, see Fig. 2.) Note that the Leapfrog scheme is
not stable for computing a diffusive evolution equation. But when using its
residue as the loss functional in N-PINN, there is no stability issue. This is
a nice feature of N-PINN which allows us to use schemes without worrying about
their conventional stability. The following subsections will discuss the
outcomes of predictions by N-PINN.

###  3.1 N-PINN with the Rusanov scheme

The first numerical scheme employed for constructing the functional form in
neural network is the Rusanov scheme which is a $1^{st}$ order conservative
scheme for conservation laws. Equation (1) can be discretized by the scheme as

$\frac{U^{n+1}_{i}-U^{n}_{i}}{\Delta
t}+\frac{\widehat{f(U)}|^{n}_{i+\frac{1}{2}}-\widehat{f(U)}|^{n}_{i-\frac{1}{2}}}{\Delta
x}=0~{},$ (7)

where

$\widehat{f(U)}|^{n}_{i+\frac{1}{2}}=\frac{1}{2}\\{{f(U^{n}_{i+1})}+{f(U^{n}_{i})}\\}+\frac{\alpha}{2}\\{U^{n}_{i}-U^{n}_{i+1}\\}~{},$
(8)

$\alpha$ is the largest characteristic speed defined by

$\alpha=\max_{i}\\{|v_{i}|+c_{i}\\}~{},$ (9)

and $c_{i}=\sqrt{\gamma\frac{p_{i}}{\rho_{i}}}$ is the sound speed. In this
N-PINN we fix $\alpha$ to a number determined by (9) using the given initial
value, multiplied by a factor, e.g., $2$. We present three examples predicted
by the N-PINN.

Example 1. Burgers’ Equation

Consider the Burgers’ equation

$\frac{\partial U}{\partial t}+\frac{\partial(\frac{1}{2}U^{2})}{\partial
x}=0,x\in[-1,1],t\in[0,1],$ (10)

with the initial condition given as

$\displaystyle U(0,x)=-\sin(\pi x)~{}.$

The viscous version of (10) is computed by neural network in Raissi et al.
[25]. The original solution with periodic boundary conditions and specified
initial condition is simulated from Chebfun package [47]. The loss function is
minimized by the L-BFGS optimizer (in TensorFlow [48, 26]) during the training
process. The total loss is defined in eqn. (4). Instead of defining the
functional form by auto-differentiation in TensorFlow, it is constructed in
(6) with residue $L$ defined to be the left-hand-side of (7).

Figure 3: Burgers’ equation: (a) to (c) exhibit the prediction and exact
solution snapshots from chosen temporal locations, $t=0.25$, $t=0.50$, and
$t=0.75$; blue solid line represents exact solution and red dashed line
represent prediction results.

Fig. 3 summarizes the prediction of solution for Burgers’ equation by N-PINN.
To proceed the training process, the initial and boundary conditions ($0$ in
this example) are applied in (5) where the initial condition is given at a set
of $256$ data points and boundary conditions are given at $200$ data points
($100$ data points for each boundary). The associated DNN consists of $10$
hidden layers and $40$ neurons in each layer. The resulting error measured by
relative $L_{2}$ norm is $1.328e-2$. Compared with the original PINN, the
error is $2$ orders of magnitude higher because Rusanov scheme is only first
order. However, the error is lower than that of the framework using the
numerical gaussian process to solve time-dependent and nonlinear partial
differential equations [49, 50].

Exampe 2. Lax Problem

The Lax problem is a Riemann problem for the Euler equations with initial
value

$\displaystyle U(0,x)=(0.445,0.311,8.928)^{T},\text{if}\ 0\leq x\leq 0.5~{},$
$\displaystyle U(0,x)=(0.5,0.0,1.4275)^{T},\text{if}\ 0.5<x\leq 1.0~{},$

where $U(t,x)=(\rho,\rho v,E)^{T}$.

Figure 4: Final time prediction $(t=0.16)$ of the Lax problem by N-PINN based
on Rusanov scheme (blue) in comparison with the reference solution (red).

Fig. 4 presents the comparative results of the Lax problem predicted by N-PINN
(based on the Rusanov scheme), and a high resolution reference solution
computed by a 3rd order finite volume scheme using non-oscillatory
hierarchical reconstruction (HR) limiting [9] and partial neighboring cells
[10] in HR, which will be used for computing reference solutions of the Euler
equations throughout the paper. Even though N-PINN is able to capture the
discontinuous solution reasonably well, it has a large numerical diffusion
around the discontinuities due to the underlying first order scheme used in
the loss function. Since the number of unknowns in the computational domain
(with $200$ cells in $[0,1]$) of the Lax problem is greater than that of the
Burgers’ equation, we increase the number of neurons accordingly. The
associated DNN consists of 6 hidden layers and 180 neurons in each layer. We
have tested multiple fully-connected neural network structures to search for
the optimal one for training, and found that the structure applied in this
case is the most robust one. The weight for the functional loss ($Err_{2}$ in
(4)) is 0.7 while the weights for initial and boundary value loss ($Err_{1}$)
are $0.15$ each. The greater weight of $Err_{2}$ is to enhance the resolution
around shocks and contacts.

Example 3. Sod Problem

The Sod problem is a Riemann problem for the Euler equations with initial
value

$\displaystyle U(0,x)=(1.0,0.0,2.5)^{T},\text{if}\ 0\leq x\leq 0.5~{},$
$\displaystyle U(0,x)=(0.125,0.0,0.25)^{T},\text{if}\ 0.5<x\leq 1.0~{}.$

Figure 5: Final time prediction $(t=0.16)$ of the Sod problem by N-PINN based
on Rusanov scheme (blue) in comparison with the reference solution (red).

Fig. 5 shows the comparison of the final-time solution of N-PINN (based on the
Rusanov scheme) and the reference solution. The numerical diffusion is large
due to the first order underlying scheme used in the loss function. The setup
of the associated DNN as well as the weights for parts of the loss function
are the same as those in the Lax problem.

###  3.2 Physics Informed Neural Network with Leapfrog Scheme and Artificial
Diffusion

Leapfrog scheme is a commonly used second order scheme for solving first order
hyperbolic PDEs. However, artificial diffusion terms of size $O(\Delta x)$
need to be added to the equations in order to damp the numerical oscillations
around discontinuities generated by a linear second order scheme. The central
discretization of the diffusion terms yields the following (unstable) scheme

$\frac{U^{n+1}_{i}-U^{n-1}_{i}}{2\Delta
t}+\frac{{f(U)}|^{n}_{i+1}-{f(U)}|^{n}_{i-1}}{2\Delta
x}-\alpha\frac{U^{n}_{i+1}-2\cdot U^{n}_{i}+U^{n}_{i-1}}{\Delta x^{2}}=0~{},$
(11)

where $\alpha=\Delta x$. Nevertheless, the residue of the scheme, i.e., the
left-hand-side of (11) can be used as residue $L$ in (6), a nice property of
the neural network. The resulted N-PINN (N-PINN based on Leapfrog and
diffusion) performs similarly to that based on the Rusanov scheme. Fig. 6
displays the final-time predicted solution of N-PINN based on Leafrog and
diffusion compared to the reference solution, which has slightly higher
resolution than that computed by N-PINN based on the Rusanov scheme. We use
the same neural network structures and weights for the loss function during
the training procedure as those for N-PINN based on the Rusanov scheme.

Figure 6: Final time prediction $(t=0.16)$ of the Lax problem by N-PINN based
on Leapfrog and diffusion (blue) in comparison with the reference solution
(red).

##  4 Numerical Physics Informed Neural Network (N-PINN) with Additional
Training Data

###  4.1 Methodology

N-PINN based on first order schemes has relatively poor resolution around
shocks and contacts because of the large numerical diffusion of the underlying
scheme. We have tried using the residue of the second order MUSCL scheme [3]
in the loss function of N-PINN but the minimization procedure has difficulty
in reducing it to the desired level. In previous studies [25, 16, 41], sparse
training data from the exact solution is used in PINN to provide more
information to the loss functional and enhance the accuracy of prediction. We
provide in the first part of the loss functional an additional (5) $21.9\%$ of
total data points ($7437$ out of $202\times 168$ space-time grid points for
the Lax problem; $3003$ out of $202\times 79$ grid points for the Sod problem)
allocated around shocks and contacts, where at the same time we decrease the
weight (e.g., from $1$ to $0.1$) of local residue error $L$ in (6).

###  4.2 Results and Discussions

We present the results of N-PINN based on Leapfrog and diffusion with
additional training data for the Lax and Sod problem. Because of the
additional training data, the number of neurons in the associated DNN
increases to $240$ each layer for the Lax problem while the number of layers
remains to be $6$. For the Sod problem, the associated DNN remains the same as
before ($6$ layers and $180$ neurons in each layer.) Meanwhile, the weights
for the loss function are adjusted. We assign the weight $0.1$ for the loss
from initial value, $0.1$ for the loss from boundary conditions, $0.5$ for the
loss from residue error of the underlying scheme, and $0.3$ for the loss from
additional training data. Fig. 7 shows the final-time profiles of density,
pressure and velocity of the N-PINN prediction of the Lax problem and
reference solution. Except for the numerical artifact around the contact for
the velocity field, the overall resolution is much improved from Fig. 6 which
has similar setup in N-PINN but without the additional training data. Fig. 8,
shows the final-time profiles of density, pressure and velocity of the N-PINN
prediction of the Sod problem and reference solution. The quality of the
prediction is much better than those in Fig. 5 which has similar setup in
N-PINN but without the additional training data.

Figure 7: Final time prediction $(t=0.16)$ of the Lax problem by N-PINN based
on Leapfrog and diffusion with additional data (blue) in comparison with the
reference solution (red).

Figure 8: Final time prediction $(t=0.16)$ of the Sod problem by N-PINN based
on Leapfrog and diffusion with additional data (blue) in comparison with the
reference solution (red).

##  5 2-Coarse-Grid Neural Network

It is an efficient strategy to use a low cost numerical solution of a PDE as
the input of a neural network to predict a more accurate solution of the PDE.
In [37], the numerical solution in a coarse mesh is used in a neural network
to predict the solution in a fine mesh. In [40], the solution of a second
order wave equation and its gradient computed by a first order scheme is used
to predict the solution computed by a higher order scheme. In this study, we
introduce a novel neural network, the 2-Coarse-Grid neural network (2CGNN) to
accurately compute solutions of systems of conservation laws which may contain
shocks and contacts. Instead of taking a spatiotemporal point $(x,t)$ as
input, 2CGNN first needs two different coarse grid solutions to be computed.
The input of the neural network is designed to incorporate a small domain of
dependence of the solution at certain space-time location as well as an
approaching sequence of approximations. Once trained, the neural network is
able to predict a fine discontinuous solution efficiently based on some low-
cost computations.

###  5.1 Input, Output and Loss Functional

Consider the scalar conservation law (1) on an interval $[a,b]$ which is
partitioned with the coarsest uniform grid $a=x_{0}<x_{1}<\dots<x_{M}=b$
having spatial grid size $\Delta x=x_{1}-x_{0}$ and time step size $\Delta t$.
Refine the grid to obtain a finer uniform grid with spatial grid size
$\frac{1}{2}\Delta x$ and time step size $\frac{1}{2}\Delta t$. Let $L$ be a
low-cost scheme used to compute (1) on both grids. Given a grid point
$x_{i^{\prime}}$ at the time $t_{n^{\prime}}$ (on the coarsest uniform grid)
where the solution is to be predicted by DNN, we choose the coarsest grid
solution (computed by $L$) at $3$ points $x_{i^{\prime}-1}$, $x_{i^{\prime}}$
and $x_{i^{\prime}+1}$ at the time level $t_{n^{\prime}-1}$ and also at the
point $(x_{i^{\prime}},t_{n^{\prime}})$ as the first part of the input, and
the finer grid solution (also computed by $L$) at the same spatiotemporal
locations as the second part of input. Note that, the chosen $4$
spatiotemporal locations of either grid enclose a (space-time) domain of
dependence of the exact solution at $(x_{i^{\prime}},t_{n^{\prime}})$ (with
$\Delta t$ satisfying the CFL restriction.) Since the two parts of the input
solution have different levels of approximation to the exact or reference
solution, DNN utilizes the information to extrapolate a prediction of the
exact or reference solution.

Denote the first part of the input as

$u^{n^{\prime}-1}_{i^{\prime}-1},u^{n^{\prime}-1}_{i^{\prime}},u^{n^{\prime}-1}_{i^{\prime}+1},u^{n^{\prime}}_{i^{\prime}}~{},$

and the second part of input as

$u^{n^{\prime\prime}-2}_{i^{\prime\prime}-2},u^{n^{\prime\prime}-2}_{i^{\prime\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}+2},u^{n^{\prime\prime}}_{i^{\prime\prime}}~{}.$

Note that the space-time index $(i^{\prime},n^{\prime})$ in the coarsest grid
refers to the same location as $(i^{\prime\prime},n^{\prime\prime})$ does in
the finer grid, $(i^{\prime}-1,n^{\prime}-1)$ refers to the same location as
$(i^{\prime\prime}-2,n^{\prime\prime}-2)$ does, and so on.

Suppose we are interested in the predicted solution at $(x,t)$ which is
referred to as $(i^{\prime},n^{\prime})$ in the coarsest grid, the input of
2CGNN is

$\\{u^{n^{\prime}-1}_{i^{\prime}-1},u^{n^{\prime}-1}_{i^{\prime}},u^{n^{\prime}-1}_{i^{\prime}+1},u^{n^{\prime}}_{i^{\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}-2},u^{n^{\prime\prime}-2}_{i^{\prime\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}+2},u^{n^{\prime\prime}}_{i^{\prime\prime}}\\}~{},$
(12)

called “input of $u$”, and the corresponding output of 2CGNN is the predicted
solution at $(x,t)$. See Fig. 9 for an illustration of 2CGNN. For the Euler
system, the input and output of 2CGNN are made up of corresponding ones for
each prime variable. For example, the input can be the vector

$\\{{\rm input}\;{\rm of}\;\rho,\;{\rm input}\;{\rm of}\;v,\;{\rm input}\;{\rm
of}\;p\\}$

with $8\times 3=24$ elements. And the corresponding output will be

$\\{\rho,v,p\\}$

with $3$ elements.

The loss function measures the difference between the output and the reference
solution corresponding to the input, and is defined as follows.

$\begin{array}[]{ccc}{\rm Loss}&=&\sum_{k}\|({\rm output}\;{\rm
corresponding}\;{\rm to}\;k^{th}\;{\rm set}\;{\rm of}\;{\rm input})-\\\
&&({\rm reference}\;{\rm solution}\;{\rm corresponding}\;{\rm
to}\;k^{th}\;{\rm set}\;{\rm of}\;{\rm input})\|_{2}^{2}~{},\end{array}$

where $\|\cdot\|_{2}$ is the $2$-norm, and the summation goes through every
set of input in the training data.

(a)

(b)

Figure 9: Schematic of 2-Coarse-Grid neural network: (a) Procedure of
transferring 2 coarse grids into input format; (b) Training procedure of
2-Coarse-Grid neural network.

###  5.2 Generation of Input and Training Data

The associated DNN for 2CGNN consists of $8$ hidden layers and each layer has
$300$ neurons. During the training process, DNN minimizes the difference
between outputs of the neural network and a reference solution by using L-BFGS
optimizer in TensorFLow (with the number of iterations of optimization
procedure under $50000$.) After the training is done, the neural network is
used to predict a solution (different from the training data) given an input
computed by the same low-cost scheme(s) and grids which are used to compute
inputs of the training data. We try to use the trained neural network to
predict final solutions for $5$ initial values of the Euler system which
include the original initial value of the Lax (or Sod) problem, $\pm 3\%$ and
$\pm 5\%$ perturbations of the initial value.

In order to generate the training data, we use a first order scheme on the
coarsest uniform grid ($50$ cells) and the finer uniform grid ($100$ cells) to
compute the input data from several initial values of the Euler system which
include $\pm 2\%$, $\pm 4\%$, $\pm 6\%$, $\pm 8\%$ and $\pm 10\%$
perturbations of the initial value of the Lax (or Sod) problem. The high
resolution reference solutions of the training data are computed on a uniform
grid with $200$ cells (note that solution values at grid points need to be
interpolated from cell averages.) The first order scheme used for generating
the inputs is a stable version of the Leapfrog and diffusion scheme (13), the
Rusanov scheme or the Leapfrog and diffusion splitting scheme (14). The time
step size for the first order schemes is fixed to a constant, e.g., $\Delta t$
for the $50$-cell grid and $\frac{1}{2}\Delta t$ for the $100$-cell grid,
bounded by the CFL restriction based on the largest characteristic speed
estimated throughout the computational domain. Using different low-cost
schemes in 2CGNN helps demonstrate the robustness of the method. Prediction
results based on schemes (13) and (14) are shown in Appendix A and B
respectively.

$\frac{U^{n+1}_{i}-U^{n-1}_{i}}{2\Delta
t}+\frac{{f(U)}|^{n}_{i+1}-{f(U)}|^{n}_{i-1}}{2\Delta
x}-\alpha\frac{U^{n-1}_{i+1}-2\cdot U^{n-1}_{i}+U^{n-1}_{i-1}}{\Delta
x^{2}}=0~{},$ (13)

where $\alpha=\Delta x$.

$\left\\{\begin{array}[]{l}\frac{\tilde{U}_{i}-U^{n-1}_{i}}{2\Delta
t}+\frac{{f(U)}|^{n}_{i+1}-{f(U)}|^{n}_{i-1}}{2\Delta x}=0~{},\\\
\frac{U^{n+1}_{i}-\tilde{U}_{i}}{\Delta
t}-\alpha\frac{\tilde{U}_{i+1}-2\cdot\tilde{U}_{i}+\tilde{U}_{i-1}}{\Delta
x^{2}}=0~{},\end{array}\right.$ (14)

where $\alpha=\Delta x$.

###  5.3 Results and Discussions for The Lax Problem

We first use the Rusanov Scheme to compute inputs of 2CGNN. The prediction of
the final-time solution of the Lax problem is shown in Fig. 10. It’s clear
that 2CGNN captures the shock and contact very sharp, and the prediction is
comparable or even better than the reference solutions used in training (and
as the “exact” solution in the figure.) Fig. 11 (12) shows the prediction of
the final-time solution of the Euler system with its initial value being
$+3\%$ ( $-3\%$) perturbation of that of the Lax problem. Fig. 13 (14) shows
the prediction of the final-time solution of the Euler system with its initial
value being $+5\%$ ( $-5\%$) perturbation of that of the Lax problem. Since
the coarsest grid for input has $50$ grid cells, the spatial grid for the
predicted solution also has $50$ grid cells. The predictions do not depict
smeared solutions like their low-cost input solutions, and are significant
improvements over those by N-PINN. Note that the cases for predictions are not
included in the training data.

Figure 10: 2CGNN prediction of final-time $(t=0.159888)$ solution of Lax
problem (dark blue), low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

Figure 11: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $+3\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

Figure 12: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $-3\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

Figure 13: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

Figure 14: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $-5\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

###  5.4 Results for The Sod Problem

The 2CGNN prediction of the final-time solution of the Sod problem is shown in
Fig. 15. It’s clear that 2CGNN captures the shock and contact really sharp,
even without the small overshoot of the reference solution (note that similar
reference solutions are used in the training of 2CGNN.) Fig. 16 shows the
prediction of the final-time solution of the Euler system with its initial
value being $+5\%$ perturbation of that of the Sod problem. The 2CGNN
predictions for the final-time solutions of the Euler system with its initial
value being $\pm 3\%$ and $-5\%$ perturbations of that of the Sod problem
perform similarly well and are not shown here. Compared to N-PINN and PINN
with additional data, 2CGNN once trained is much more efficient and accurate
to predict the solution of the Euler system as long as it doesn’t vary too
much from the training data.

Figure 15: 2CGNN prediction of final-time $(t=0.158213)$ solution (dark blue)
of the Sod problem, low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

Figure 16: 2CGNN prediction of final-time $(t=0.158213)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Sod problem, low-fidelity input solutions (blue and red) by Rusanov
scheme on $2$ different grids (with $50$ and $100$ cells resp.), and “exact”
(reference) solution (green).

###  5.5 Variations of the Input

There could be alternative input formats. For example, one could use as many
data points as possible from the finer grid, and the input for a scalar
conservation law can be changed from (12) to

$\\{u^{n^{\prime}-1}_{i^{\prime}-1},u^{n^{\prime}-1}_{i^{\prime}},u^{n^{\prime}-1}_{i^{\prime}+1},u^{n^{\prime}}_{i^{\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}-2},u^{n^{\prime\prime}-2}_{i^{\prime\prime}-1},u^{n^{\prime\prime}-2}_{i^{\prime\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}+1},u^{n^{\prime\prime}-2}_{i^{\prime\prime}+2},u^{n^{\prime\prime}}_{i^{\prime\prime}}\\}~{},$
(15)

similarly for the Euler system. Prediction results from this type of input is
similar to (slightly better than) those from input (12). Another variation is
to change input (15) to

$\\{u^{n^{\prime}-1}_{i^{\prime}-1},u^{n^{\prime}-1}_{i^{\prime}},u^{n^{\prime}-1}_{i^{\prime}+1},u^{n^{\prime}}_{i^{\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}-1},u^{n^{\prime\prime}-2}_{i^{\prime\prime}},u^{n^{\prime\prime}-2}_{i^{\prime\prime}+1},u^{n^{\prime\prime}}_{i^{\prime\prime}}\\}~{},$
(16)

similarly for the Euler system. Prediction results from this type of input is
similar to (slightly worse than) those from input (12).

Input (16) computed by the Leapfrog diffusion and splitting scheme (14)
generates some artifacts for predictions of the Sod problem, see Fig. 17.

Figure 17: 2CGNN (with non-standard input (16)) prediction of final-time
$(t=0.158213)$ solution (dark blue) of the Sod problem, low-fidelity input
solutions (blue and red) by Leapfrog and diffusion splitting scheme (14) on
$2$ different grids (with $50$ and $100$ cells resp.), and “exact” (reference)
solution (green).

We use $9$ layers and $100$ neurons on each layer for 2CGNN (with non-standard
input (15)) to predict the Sod problem while there is no change on neurons and
layers for other cases. Table 1 summarizes the prediction results of input
(15) computed by the leapfrog diffusion splitting scheme (14).

Table 1: Relative $l_{2}$ errors of 2CGNN (with non-standard input (15)) predictions of solutions on coarsest uniform spatial grid (50 cells) of the Euler system (with initial values perturbed from those of the Lax or Sod problem) with low-fidelity input solutions computed by the Leapfrog and diffusion splitting scheme (14) on $2$ different grids (with $50$ and $100$ cells.) | Lax Problem | Sod Problem
---|---|---
Initial Value | Relative $l_{2}$ Error | Relative $l_{2}$ Error
Original | $8.17e-03$ | $8.35e-03$
$+3\%$ | $4.06e-03$ | $2.81e-02$
$-3\%$ | $5.96e-03$ | $6.63e-03$
$+5\%$ | $5.27e-03$ | $9.37e-03$
$-5\%$ | $1.06e-02$ | $9.47e-03$
$+7\%$ | $8.51e-03$ | $8.27e-03$
$-7\%$ | $1.49e-02$ | $2.46e-02$

##  6 2-Diffusion-Coefficient Neural Network

Instead of computing for the input on two grids, one could generate the input
on a single grid but for two perturbed equations with two different diffusion
coefficients. This is referred to as the 2-Diffusion-Coefficient Neural
Network (2DCNN.) The perturbed scalar conservation law can be written as

$\frac{\partial U}{\partial t}+\frac{\partial f(U)}{\partial
x}=\alpha\frac{\partial^{2}U}{\partial
x^{2}},\;x\in\Omega\subset\mathcal{R},\;t\in[0,T]~{}.$ (17)

To fill the two parts of an input, one could compute equation (17) with two
different values of $\alpha$ using a low-cost scheme on a single grid. For
example, we can approximate (17) using the Leapfrog and diffusion splitting
scheme (14) with $\alpha=\Delta x$ and $c\Delta x$ as follows

$\left\\{\begin{array}[]{l}\frac{\tilde{U}_{i}-U^{n-1}_{i}}{2\Delta
t}+\frac{{f(U)}|^{n}_{i+1}-{f(U)}|^{n}_{i-1}}{2\Delta x}=0~{},\\\
\frac{U^{n+1}_{i}-\tilde{U}_{i}}{\Delta t}-\Delta
x\frac{\tilde{U}_{i+1}-2\cdot\tilde{U}_{i}+\tilde{U}_{i-1}}{\Delta
x^{2}}=0~{},\end{array}\right.$ (18)

and

$\left\\{\begin{array}[]{l}\frac{\tilde{V}_{i}-V^{n-1}_{i}}{2\Delta
t}+\frac{{f(V)}|^{n}_{i+1}-{f(V)}|^{n}_{i-1}}{2\Delta x}=0~{},\\\
\frac{V^{n+1}_{i}-\tilde{V}_{i}}{\Delta t}-c\Delta
x\frac{\tilde{V}_{i+1}-2\cdot\tilde{V}_{i}+\tilde{V}_{i-1}}{\Delta
x^{2}}=0~{}.\end{array}\right.$ (19)

Then the input computed on a single uniform grid can be written as

$\\{U^{n-1}_{i-1},U^{n-1}_{i},U^{n-1}_{i+1},U^{n}_{i},V^{n-1}_{i-1},V^{n-1}_{i},V^{n-1}_{i+1},V^{n}_{i}\\}~{}.$
(20)

similarly for the Euler system. See Fig. 18 and 19 for the predictions by
2DCNN, with the input computed on a uniform grid with $100$ cells and $c=4$.
Note that for the second equation of (19), the larger diffusion coefficient
may require a time step size smaller than that of the first equation of (19).
Therefore one may need to compute the second equation of (19) in two time
steps, with the time step size $\frac{1}{2}\Delta t$ for each time step.

Figure 18: 2DCNN prediction of final-time $(t=0.159888)$ solution of Lax
problem (dark blue), low-fidelity input solutions (blue and red) by the
Leapfrog and diffusion splitting scheme (14) on the uniform grid with $100$
cells, with diffusion coefficient $\alpha=\Delta x$ and $4\Delta x$ resp., and
“exact” (reference) solution (green).

Figure 19: 2DCNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by the
Leapfrog and diffusion splitting scheme (14) on the uniform grid with $100$
cells, with diffusion coefficient $\alpha=\Delta x$ and $4\Delta x$ resp., and
“exact” (reference) solution (green).

##  7 Summary and Comparison of All Methods

The tables below summarize the results of all methods we have tested in
aforementioned sections.

Table 2: Comparison of N-PINN: Relative $l_{2}$ norm errors between the predicted and “exact” (reference) solutions for the Lax problem, Sod problem and Burgers’ equation. N-PINN with | Lax Problem | Sod Problem | Burgers’
---|---|---|---
|  |  | Equation
| Density | Pressure | Velocity | Density | Pressure | Velocity |
Rusanov Scheme | $1.43e-01$ | $3.85e-02$ | $8.63e-02$ | $6.06e-02$ | $5.63e-02$ | $2.76e-01$ | $1.33e-02$
Leapfrog and Diffusion | $1.67e-01$ | $4.63e-02$ | $9.07e-02$ | $4.75e-02$ | $4.94e-02$ | $2.25e-01$ | $1.11e-02$
Leapfrog and Diffusion |  |  |  |  |  |  |
w/ Additional Data | $4.33e-02$ | $7.94e-03$ | $2.79e-02$ | $1.26e-02$ | $1.56e-02$ | $1.81e-02$ | N/A
Table 3: Relative $l_{2}$ errors of 2CGNN predictions of solutions on coarsest uniform spatial grid (50 cells) of the Euler system (with initial values perturbed from those of the Lax or Sod problem) with low-fidelity input solutions computed by the Rusanov scheme on $2$ different grids (with $50$ and $100$ cells.) | Lax Problem | Sod Problem
---|---|---
Initial Value | Relative $l_{2}$ Error | Relative $l_{2}$ Error
Original | $5.97e-03$ | $6.20e-03$
$+3\%$ | $7.29e-03$ | $5.66e-03$
$-3\%$ | $9.10e-03$ | $9.05e-03$
$+5\%$ | $8.61e-03$ | $7.84e-03$
$-5\%$ | $1.32e-02$ | $1.29e-02$
$+7\%$ | $9.96e-03$ | $9.72e-03$
$-7\%$ | $1.75e-02$ | $1.45e-02$
Table 4: Relative $l_{2}$ errors of 2CGNN predictions of solutions on coarsest uniform spatial grid (50 cells) of the Euler system (with initial values perturbed from those of the Lax or Sod problem) with low-fidelity input solutions computed by the stabilized Leapfrog and diffusion scheme (13) on $2$ different grids (with $50$ and $100$ cells.) | Lax Problem | Sod Problem
---|---|---
Initial Value | Relative $l_{2}$ Error | Relative $l_{2}$ Error
Original | $8.73e-03$ | $6.22e-03$
$+3\%$ | $5.47e-03$ | $7.68e-03$
$-3\%$ | $1.14e-02$ | $1.13e-02$
$+5\%$ | $1.01e-02$ | $6.47e-03$
$-5\%$ | $1.29e-02$ | $1.11e-02$
$+7\%$ | $1.42e-02$ | $8.37e-03$
$-7\%$ | $1.53e-02$ | $1.77e-02$
Table 5: Relative $l_{2}$ errors of 2CGNN predictions of solutions on coarsest uniform spatial grid (50 cells) of the Euler system (with initial values perturbed from those of the Lax or Sod problem) with low-fidelity input solutions computed by the Leapfrog and diffusion splitting scheme (14) on $2$ different grids (with $50$ and $100$ cells.)  | Lax Problem | Sod Problem
---|---|---
Initial Value | Relative $l_{2}$ Error | Relative $l_{2}$ Error
Original | $8.56e-03$ | $5.72e-03$
$+3\%$ | $4.80e-03$ | $8.03e-03$
$-3\%$ | $6.86e-03$ | $8.26e-03$
$+5\%$ | $4.96e-03$ | $6.60e-03$
$-5\%$ | $7.05e-03$ | $1.84e-02$
$+7\%$ | $6.22e-03$ | $8.55e-03$
$-7\%$ | $8.28e-03$ | $1.05e-02$
Table 6: Relative $l_{2}$ errors of 2DCNN predictions of solutions on finer uniform spatial grid (100 cells) of the Euler system (with initial values perturbed from those of the Lax or Sod problem) with low-fidelity input solutions computed by the Leapfrog and diffusion splitting scheme (14) with $2$ different diffusion coefficients ($\alpha=\Delta x$ and $4\Delta x$) on a uniform grid with $100$ cells.  | Lax Problem | Sod Problem
---|---|---
Initial Value | Relative $l_{2}$ Error | Relative $l_{2}$ Error
Original | $5.24e-03$ | $6.47e-03$
$+3\%$ | $3.43e-03$ | $6.80e-03$
$-3\%$ | $5.34e-03$ | $8.57e-03$
$+5\%$ | $7.13e-03$ | $7.52e-03$
$-5\%$ | $6.59e-03$ | $1.07e-02$
$+7\%$ | $8.19e-03$ | $8.87e-03$
$-7\%$ | $7.51e-03$ | $1.30e-02$

##  8 Conclusion

We have studied a few neural networks based on first order schemes to solve
the Burgers’ equation and one dimensional Euler system. An N-PINN is able to
generate an approximate solution comparable to those obtained from first order
schemes. Given additional data of the solution, N-PINN can improve its
approximation conveniently. We also introduce 2CGNN and 2DCNN which
significantly improve from N-PINN and are efficient to use after they have
been trained. In fact, they are able to output very sharp shocks and contacts
out of smeared solution profiles in the input. And they don’t seem to be very
sensitive to the low-cost schemes used for computing inputs either, see
Appendix A and B for additional prediction results with some alternative
schemes used for computing inputs. We will extend 2CGNN and 2DCNN to
multidimensional nonlinear hyperbolic systems in the future.

## Acknowledgements

The authors thank Prof. Wen Shen, Penn. State U. for a helpful discussion.

## Appendix A. Stabilized Leapfrog and Diffusion Scheme (13) for Computing
Inputs for 2CGNN

Figure 20: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Lax problem, low-fidelity input solutions (blue and red) by stabilized
Leapfrog and diffusion scheme (13) on $2$ different grids (with $50$ and $100$
cells resp.), and “exact” (reference) solution (green).

Figure 21: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by stabilized
Leapfrog and diffusion scheme (13) on $2$ different grids (with $50$ and $100$
cells resp.), and “exact” (reference) solution (green).

Figure 22: 2CGNN prediction of final-time $(t=0.158213)$ solution (dark blue)
of the Sod problem, low-fidelity input solutions (blue and red) by stabilized
Leapfrog and diffusion scheme (13) on $2$ different grids (with $50$ and $100$
cells resp.), and “exact” (reference) solution (green).

Figure 23: 2CGNN prediction of final-time $(t=0.158213)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Sod problem, low-fidelity input solutions (blue and red) by stabilized
Leapfrog and diffusion scheme (13) on $2$ different grids (with $50$ and $100$
cells resp.), and “exact” (reference) solution (green).

## Appendix B. Leapfrog and Diffusion Splitting Scheme (14) for Computing
Inputs for 2CGNN

Figure 24: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Lax problem, low-fidelity input solutions (blue and red) by Leapfrog
and diffusion splitting scheme (14) on $2$ different grids (with $50$ and
$100$ cells resp.), and “exact” (reference) solution (green).

Figure 25: 2CGNN prediction of final-time $(t=0.159888)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Lax problem, low-fidelity input solutions (blue and red) by Leapfrog
and diffusion splitting scheme (14) on $2$ different grids (with $50$ and
$100$ cells resp.), and “exact” (reference) solution (green).

Figure 26: 2CGNN prediction of final-time $(t=0.158213)$ solution (dark blue)
of the Sod problem, low-fidelity input solutions (blue and red) by Leapfrog
and diffusion splitting scheme (14) on $2$ different grids (with $50$ and
$100$ cells resp.), and “exact” (reference) solution (green).

Figure 27: 2CGNN prediction of final-time $(t=0.158213)$ solution (dark blue)
of the Euler system with its initial value being $+5\%$ perturbation of that
of the Sod problem, low-fidelity input solutions (blue and red) by Leapfrog
and diffusion splitting scheme (14) on $2$ different grids (with $50$ and
$100$ cells resp.), and “exact” (reference) solution (green).

## References

  * [1] S. K. Godunov, A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations, Mat. Sb. 47 (1959) 271–306.
  * [2] B. van Leer, Towards the ultimate conservative difference scheme I, Lecture Notes in Phys. 18 (1973) 163–168.
  * [3] B. van Leer, Towards the ultimate conservative difference scheme v, a second-order sequel to godunov’s method, J. Comput. Phys. 32 (1979) 101–136.
  * [4] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order accuracy essentially non-oscillatory schemes iii, J. Comput. Phys. 71 (2) (1987) 231–303.
  * [5] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471.
  * [6] X. D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1) (1994) 200–212.
  * [7] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228.
  * [8] Y. Liu, C.-W. Shu, E. Tadmor, M.-P. Zhang, Central discontinuous Galerkin methods on overlapping cells with a non-oscillatory hierarchical reconstruction, SIAM J. Num. Anal. 45 (2007) 2442–2467.
  * [9] Y. Liu, C.-W. Shu, E. Tadmor, M.-P. Zhang, Non-oscillatory hierarchical reconstruction for central and finite volume schemes, Communications in Computational Physics 2 (2007) 933–963.
  * [10] Z. Xu, Y. Liu, C.-W. Shu, Hierarchical reconstruction for discontinuous Galerkin methods on unstructured grids with a weno type linear reconstruction and partial neighboring cells, J. Comput. Phys. 228 (2009) 2194–2212.
  * [11] V. Yang, Modeling of supercritical vaporization, mixing, and combustion processes in liquid-fueled propulsion systems, Proceedings of the Combustion Institute 28 (1) (2000) 925–942.
  * [12] X. Wang, V. Yang, Supercritical mixing and combustion of liquid-oxygen/ kerosene bi-swirl injectors, Journal of Propulsion and Power 33 (2) (2017) 316–322.
  * [13] U. Unnikrishnan, H. Huo, X. Wang, V. Yang, Subgrid scale modeling considerations for large eddy simulation of supercritical turbulent mixing and combustion, Physics of Fluids 33 (7) (2021) 075112.
  * [14] R. Teyssier, B. Commerçon, Numerical methods for simulating star formation, Frontiers in Astronomy and Space Sciences 6 (2019) 51.
  * [15] A. Darcy, A. Louie, L. Roberts, Machine learning and the profession of medicine, JAMA 315 (6) (2016) 551–552.
  * [16] G. Karniadakis, I. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (2021) 422–440.
  * [17] B. Peherstorfer, K. Willcox, M. Gunzburger, Survey of multifidelity methods in uncertainty propagation, inference, and optimization, SIAM Review 60 (3) (2018) 550–591.
  * [18] A. Krizhevsky, I. Sutskever, G. Hinton, Imagenet classification with deep convolutional neural networks, Vol. 2, 2012, pp. 1097–1105, cited By 49543.
  * [19] Y. Chang, X. Wang, L. Zhang, Y. Li, S. Mak, C. Wu, V. Yang, Common kernel-smoothed proper orthogonal decomposition (ckspod): An efficient reduced-order model for emulation of spatiotemporally evolving flow dynamics (2021).  arXiv:2101.08893.
  * [20] Y. Chang, L. Zhang, X. Wang, S. Yeh, S. Mak, C. Sung, C. Wu, V. Yang, Kernel-smoothed proper orthogonal decomposition–based emulation for spatiotemporally evolving flow dynamics prediction, AIAA Journal 57 (12) (2019) 5269–5280.
  * [21] S. Mak, C. Sung, X. Wang, S. Yeh, Y. Chang, V. Joseph, V. Yang, C. F. J. Wu, An efficient surrogate model for emulation and physics extraction of large eddy simulations, Journal of the American Statistical Association 113 (524) (2018) 1443–1456.
  * [22] S. Xin, S. Nousias, K. N. Kutulakos, A. C. Sankaranarayanan, S. G. Narasimhan, I. Gkioulekas, A theory of fermat paths for non-line-of-sight shape reconstruction, in: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2019, pp. 6793–6802.
  * [23] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, I. Polosukhin, Attention is all you need, 2017.
  * [24] L. Chen, R. Paleja, M. Gombolay, Learning from suboptimal demonstration via self-supervised reward regression (2020).  arXiv:2010.11723.
  * [25] M. Raissi, P. Perdikaris, G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.
  * [26] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: A survey, J. Mach. Learn. Res. 18 (1) (2017) 5595–5637.
  * [27] X. Jin, S. Cai, H. Li, G. Karniadakis, Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations, Journal of Computational Physics 426 (2021) 109951.
  * [28] P. C. D. Leoni, T. Zaki, G. Karniadakis, C. Meneveau, Two-point stress–strain-rate correlation structure and non-local eddy viscosity in turbulent flows, Journal of Fluid Mechanics 914 (2021) A6.
  * [29] M. Raissi, H. Babaee, P. Givi, Deep learning of turbulent scalar mixing, Phys. Rev. Fluids 4 (2019) 124501.
  * [30] M. Raissi, Z. Wang, M. Triantafyllou, G. Karniadakis, Deep learning of vortex-induced vibrations, Journal of Fluid Mechanics 861 (2019) 119–137.
  * [31] M. Raissi, A. Yazdani, G. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367 (6481) (2020) 1026–1030.
  * [32] Q. Lou, X. Meng, G. Karniadakis, Physics-informed neural networks for solving forward and inverse flow problems via the boltzmann-bgk formulation (2020).  arXiv:2010.09147.
  * [33] Y. Chen, L. Lu, G. Karniadakis, L. D. Negro, Physics-informed neural networks for inverse problems in nano-optics and metamaterials, Opt. Express 28 (8) (2020) 11618–11633.
  * [34] M. Mahmoudabadbozchelou, M. Caggioni, S. Shahsavari, W. Hartt, G. Karniadakis, S. Jamali, Data-driven physics-informed constitutive metamodeling of complex fluids: A multifidelity neural network (mfnn) framework, Journal of Rheology 65 (2) (2021) 179–198.
  * [35] S. Cai, Z. Wang, S. Wang, P. Perdikaris, G. Karniadakis, Physics-Informed Neural Networks for Heat Transfer Problems, Journal of Heat Transfer 143 (6), 060801.
  * [36] Z. Mao, L. Lu, O. Marxen, T. A. Zaki, G. Karniadakis, Deepm& mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators (2020).  arXiv:2011.03349.
  * [37] D. Liu, Y. Wang, Multi-fidelity physics-constrained neural network and its application in materials modeling, Journal of Mechanical Design 141 (12), 121403\.
  * [38] H. Gao, L. Sun, J. Wang, Phygeonet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state pdes on irregular domain, Journal of Computational Physics 428 (2021) 110079.
  * [39] A. F. Psaros, K. Kawaguchi, G. Karniadakis, Meta-learning pinn loss functions (2021).  arXiv:2107.05544.
  * [40] H. Nguyen, R. Tsai, A stable parareal-like method for the second order wave equation, J. Comput. Phys. 405 (2020) 109156.
  * [41] Z. Mao, A. D. Jagtap, G. Karniadakis, Physics-informed neural networks for high-speed flows, Computer Methods in Applied Mechanics and Engineering 360 (2020) 112789.
  * [42] C. Michoski, M. Milosavljević, T. Oliver, D. R. Hatch, Solving differential equations using deep neural networks, Neurocomputing 399 (2020) 193–212.
  * [43] A. Harten, Eno schemes with subcell resolution, J. Comput. Phys. 83 (1) (1989) 48–184.
  * [44] W. Shen, M. R. Park, Optimal tracing of viscous shocks in solutions ofviscous conservation laws, SIAM J. MATH. ANAL. 38 (5) (2007) 1474–1488.
  * [45] K. Hornik, M. Stinchcombe, H. White, Multilayer feedforward networks are universal approximators, Neural Networks 2 (5) (1989) 359–366.
  * [46] V. V. Rusanov, Calculation of intersection of non-steady shock waves with obstacles, J. Comput. Math. Phys. USSR 1 (1961) 267–279.
  * [47] T. A. Driscoll, N. Hale, L. N. Trefethen, Chebfun Guide, Pafnuty Publications, 2014\.
  * [48] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. D. et al., Tensorflow: Large-scale machine learning on heterogeneous distributed systems (2016).  arXiv:1603.04467.
  * [49] M. Raissi, P. Perdikaris, G. Karniadakis, Inferring solutions of differential equations using noisy multi-fidelity data, Journal of Computational Physics 335 (2017) 736–746.
  * [50] M. Raissi, P. Perdikaris, G. Karniadakis, Numerical gaussian processes for time-dependent and nonlinear partial differential equations, SIAM Journal on Scientific Computing 40 (1) (2018) A172–A198.