Energy-Dissipative Evolutionary Deep Operator Neural Networks
Energy-Dissipative Evolutionary Deep Operator Neural Network is an operator
learning neural network. It is designed to seed numerical solutions for a class
of partial differential equations instead of a single partial differential
equation, such as partial differential equations with different parameters or
different initial conditions. The network consists of two sub-networks, the
Branch net and the Trunk net. For an objective operator G, the Branch net
encodes different input functions u at the same number of sensors, and the
Trunk net evaluates the output function at any location. By minimizing the
error between the evaluated output q and the expected output G(u)(y), DeepONet
generates a good approximation of the operator G. In order to preserve
essential physical properties of PDEs, such as the Energy Dissipation Law, we
adopt a scalar auxiliary variable approach to generate the minimization
problem. It introduces a modified energy and enables unconditional energy
dissipation law at the discrete level. By taking the parameter as a function of
time t, this network can predict the accurate solution at any further time with
feeding data only at the initial state. The data needed can be generated by the
initial conditions, which are readily available. In order to validate the
accuracy and efficiency of our neural networks, we provide numerical
simulations of several partial differential equations, including heat
equations, parametric heat equations and Allen-Cahn equations.
Jiahao Zhang, Shiheng Zhang, Jie Shen, Guang Lin
# Energy-Dissipative Evolutionary Deep Operator Neural Networks

###### Abstract

Energy-Dissipative Evolutionary Deep Operator Neural Network is an operator learning neural network. It is designed to seek numerical solutions for a class of partial differential equations instead of a single partial differential equation, such as partial differential equations with different parameters or different initial conditions. The network consists of two sub-networks, the Branch net, and the Trunk net. For an objective operator \(\mathcal{G}\), the Branch net encodes different input functions \(u\) at the same number of sensors \(y_{i},i=1,2,\cdots,m\), and the Trunk net evaluates the output function at any location. By minimizing the error between the evaluated output \(q\) and the expected output \(\mathcal{G}(u)(y)\), DeepONet generates a good approximation of the operator \(\mathcal{G}\). In order to preserve essential physical properties of PDEs, such as the Energy Dissipation Law, we adopt a scalar auxiliary variable approach to generate the minimization problem. It introduces a modified energy and enables unconditional energy dissipation law in the discrete level. By taking the parameter as a function of the time \(t\) variable, this network can predict the accurate solution at any further time with feeding data only at the initial state. The data needed can be generated by the initial conditions, which are readily available. In order to validate the accuracy and efficiency of our neural networks, we provide numerical simulations of several partial differential equations, including heat equations, parametric heat equations, and Allen-Cahn equations.

A Article history:

Operator Learning Evolutionary Neural Networks Energy Dissipative Parametric equation Scalar auxiliary variable Deep learning

## 1 Introduction

Operator learning is a popular and challenging problem with potential applications across various disciplines. The opportunity to learn an operator over a domain in Euclidean spaces[1] and Banach spaces[2] opens a new class of problems in neural network design with generalized applicability. In application to solve partial differential equations(PDEs), operator learning has the potential to predict accurate solutions for the PDE by acquiring extensive prior
knowledge [3; 4; 5; 6; 7; 8; 9; 10; 11]. In a recent paper[12], Lu, Jin, and Karniadakis proposed an operator learning method with some deep operator networks, named as DeepONets. It is based on the universal approximation theorem [13; 14; 15]. The goal of this neural network is to learn an operator instead of a single function, which is usually the solution of a PDE. For any operator \(\mathcal{G}\) on a domain \(\Omega\), we can define \(\mathcal{G}\) as a mapping from \(\Omega^{*}\rightarrow\Omega^{*}\) with \(\mathcal{G}(u)(y)\in R\) for any \(y\in\Omega\). \(\mathcal{G}(u)(y)\) is the expected output of the neural network, which is usually a real number. The objective of the training is to obtain an approximation of \(\mathcal{G}\), where we need to represent operators and functions in a discrete form. In practice, it is very common to represent a continuous function or operator by the values evaluated at finite and enough locations \(\{x_{1},x_{2},\cdots,x_{m}\}\), which is called "sensors" in DeepONet. The network takes \([u(x_{1}),u(x_{2}),\cdots,u(x_{m})]\) and \(y\) as the input. The loss function is the difference between the output \(q\) and the expected output \(\mathcal{G}(u)(y)\). Generally, there are two kinds of DeepONet, Stacked DeepONet, and Unstacked DeepONet. The Stacked DeepONet consists of \(p\) branch networks and one trunk network. The number of the Trunk networks of the Unstacked DeepONet is the same as the DeepONet, but the Unstacked DeepONet merges all the \(p\) branch networks into a single one. An Unstacked DeepONet combines two sub-networks, Branch net, and Trunk net. The Branch net encodes the input function \(u\) at some sensors, \(\{x_{i}\in\Omega\,|\,i=1,\cdots,m\}\). The output of the Branch net consists of \(p\) neurons, where each neuron can be seen as a scalar, \(b_{j}=b_{j}(u(x_{1}),u(x_{2}),\cdots,u(x_{m}))\), \(j=1,2,\cdots,p\). The Trunk net encodes some evaluation points \(\{y_{k}\in\Omega|k=1,\cdots,n\}\), while the output also consists of \(p\) neurons and each neuron is a scalar \(g_{j}=g_{j}(y_{1},y_{2},\cdots,y_{n})\), \(j=1,2,\cdots,p\). The evaluation point \(y_{i}\) can be arbitrary in order to obtain the loss function. The number of neurons at the last layer of the Trunk net and the Branch net is the same. Hence, the output of the DeepONet can be written as an inner product of \((b_{1},b_{2},\cdots,b_{p})\) and \((g_{1},g_{2},\cdots,g_{p})\). In other words, the relationship between the expected output and the evaluated output is \(\mathcal{G}(u)(y)\approx\sum_{j=1}^{p}b_{j}g_{j}\). The DeepONet is an application of the Universal Approximation Theorem for Operator, which is proposed by Chen & Chen [16]:

**Theorem 1.1** (Universal Approximation Theorem for Operator).: _Suppose that \(\Omega_{1}\) is a compact set in \(X\), \(X\) is a Banach Space, \(V\) is a compact set in \(C(\Omega_{1})\), \(\Omega_{2}\) is a compact set in \(\boldsymbol{R}^{d}\), \(\sigma\) is a continuous non-polynomial function, \(\mathcal{G}\) is a nonlinear continuous operator, which maps \(v\) into \(C(\Omega_{2})\), then for any \(\epsilon>0\), there are positive integers \(M,N,m\), constants \(c_{i}^{k},\zeta_{k},\xi_{ij}^{k}\in\boldsymbol{R}\), points \(\omega_{k}\in\boldsymbol{R}^{n},x_{j}\in K_{1},i=1,\cdots,M\), \(k=1,\cdots,N,j=1,\cdots,m\), such that_

\[|\,\mathcal{G}(u)(y)-\sum_{k=1}^{N}\sum_{i=1}^{M}c_{i}^{k}\sigma\left(\sum_{j =1}^{m}\xi_{ij}^{k}u\left(x_{j}\right)+\theta_{i}^{k}\right)\cdot\sigma\left( \omega_{k}\cdot y+\zeta_{k}\right)|<\epsilon\]

_holds for all \(u\in V\) and \(y\in\Omega_{2}\)._

For any time-dependent PDE, the training data is the form of \((u,y,\mathcal{G}(u)(y))\), where \(u\) in the discrete form can be represented as \([u(x_{1}),u(x_{2}),\cdots,u(x_{m})]\) in the neural network. In the original paper, they used the classic FNN[17] as the baseline model. For dynamic systems, various network architectures are used, including residual networks[18], convolutional NNs(CNNs)[19; 20], recurrent NNs(RNNs)[21], neural jump stochastic differential equations[22] and neural ordinary differential equations[23]. The training performance is very promising. It predicts accurate solutions of many nonlinear ODEs and PDEs, including the simple dynamic system, gravity pendulum system, and diffusion-reaction system. However, the training data need to be generated at each time step, so it is very expensive to train the network. For a lot of initial value problems, there is no any information of \(u(x,t)\) except \(t=0\). It is very natural to raise a question: _Can we learn an operator of a kind of time-dependent PDEs with only initial conditions?_

Inspired by the Evolutionary Deep Neural Network(EDNN)[24], it is more convenient to learn an operator at a fixed time instead of an operator with not only spatial variables but also a time variable. With the loss of generality, we can take the time variable \(t\) to be \(0\) in initial value problems. Once obtained the operator at the initial time, many traditional numerical methods can be used to update the solution. More specifically, assuming that the initial condition operator has been trained well, we can consider the parameters of the Branch net and the Trunk net as a function with respect to the time variable as shown in Figure 1. More specifically, for a given initial value problem,

\[\begin{cases}\dfrac{\partial u}{\partial t}=s(u)\\ u(x,0)=f(x),\quad x\in\omega\end{cases} \tag{1}\]

the objective is to approximate the operator \(\mathcal{G}:u\mapsto\mathcal{G}(u)\). The input is \(([u(x_{1}),u(x_{2}),\cdots,u(x_{m})],y,\mathcal{G}(u)(y))\), where \([x_{1},x_{2},\cdots,x_{m}]\) are the sensors and \(\mathcal{G}(u)(y)=f(y)\). The training process at the initial step is the same as the DeepONet,
so we can use the same architecture to train the initial condition operator. The output of the Branch net can be written as \(\mathbf{b}=\mathbf{b}(u(x_{1},0),u(x_{2},0),\cdots,u(x_{m_{1}},0))=\mathbf{b}^{\prime}(x_{1}, x_{2},\cdots,x_{m};W_{1})\), where \(W_{1}\) are the parameters in the Branch net. The output of the Trunk net can be written as \(\mathbf{g}=\mathbf{g}(y;W_{2})\), where \(W_{2}\) are the parameters in the Trunk net. Once trained well, we will regard the parameters as a function of \(t\) and \(W_{1}\), \(W_{2}\) as the initial conditions of \(W_{1}(t)\) and \(W_{2}(t)\). By the architecture of the Unstacked DeepONet, we can write the solution at initial time \(t_{0}=0\) as

\[u(x,t_{0})\approx\sum_{j=1}^{p}b_{j}g_{j}=\mathbf{b}^{T}\mathbf{g}\text{ for any given initial condition }f(x) \tag{2}\]

We do not need any more data to obtain the approximation of \(u(x,t_{1})\). \(u(x,t_{1})\) should be consistent with \(W_{1}(t_{1})\) and \(W_{2}(t_{1})\). With the idea of the numerical solver for PDEs, it is easy to obtain \(W_{1}(t_{1})\) and \(W_{2}(t_{1})\) if \(\frac{\partial W_{1}}{\partial t}\) and \(\frac{\partial W_{2}}{\partial t}\) are known. The time derivative of the solution \(u\) can be written by a chain rule:

\[\frac{\partial u}{\partial t}=\frac{\partial u}{\partial W}\frac{\partial W}{ \partial t} \tag{3}\]

where \(W\) consists of \(W_{1}\) and \(W_{2}\). \(\frac{\partial W}{\partial t}\) can be solved by a least square problem. Once we get \(\frac{\partial W}{\partial t}\), we can use any traditional time discretization schemes to get \(W^{n+1}\) with \(W^{n}\).

Fig. 1: Energy-Dissipative Evolutionary Deep Operator Neural Network. The yellow block represents input at sensors and the blue block represents subnetworks. The green blocks represent the output of the subnetworks and also the last layer of the EDE-DeepONet. The difference between the stacked and unstacked EDE-DeepONet is the number of Branch nets. In the right minimization problem, the energy term \(r^{2}\) can be shown to be dissipative, i.e. \((r^{n+1})^{2}\leq(r^{n})^{2}\), where \(\mathcal{J}(\gamma_{1},\gamma_{2})=\frac{1}{2}\left\|\sum_{k=1}^{p}\frac{ \partial g_{k}(W_{1}^{2})}{\partial W_{1}^{2}}\gamma_{1}\mathbf{b}_{k}(W_{2}^{n} )+\sum_{k=1}^{p}g_{k}(W_{1}^{n})\frac{\partial g_{k}(W_{2}^{n})}{\partial W_{ 2}^{n}}\gamma_{2}-\frac{r^{n+1}}{\sqrt{E(w^{n})}}N_{x}(\mathbf{u}^{n})\right\|_{2}^ {2}\).
The choice of the traditional time discretization scheme is dependent on the specific problem. The Euler or Runge-Kutta methods are commonly used in the evolutionary network. We are going to introduce a method with unconditional energy dissipation, which is the Energy-Dissipative Evolutionary Deep Operator Neural Network(EDE-DeepONet). Many kinds of PDEs are derived from basic physical laws, such as Netwon's Law, Conservation Law and Energy Dissipation Law. In many areas of science and engineering, particularly in the field of materials science, gradient flows are commonly employed in mathematical models[25; 26; 27; 28; 29; 30; 31; 32]. When approximating the solution of a certain PDE, it is desirable to satisfy these laws. We consider a gradient flow problem,

\[\frac{\partial u}{\partial t}=-\frac{\delta E}{\delta u}, \tag{4}\]

where \(E\) is a certain free energy functional. Since the general explicit Euler method does not possess the unconditionally dissipative energy dissipation law, we applied a scalar auxiliary variable(SAV) method[33] to generate the required least square problem. It introduces a new modified energy and the unconditionally dissipative modified energy dissipation law is satisfied for each iterative step. SAV method has been applied to solve plenty of PDEs with thermodynamically consistent property. It is robust, easy to implement and accurate to predict the solution. Introducing this method to neural network helps us explore how to combine neural network models and physical laws.

The objectives of this article is:

* Designing an operator learning neural network without data except the given information.
* Predicting solutions of parametric PDEs after a long time period.
* Keeping energy dissipative property of a dynamic system.

Our main contributions are:

* Constructing an evolutionary operator learning neural network to solve PDEs.
* Solving a kind of PDEs with different parameters in a single neural network.
* Introducing the modified energy in the neural network and applying SAV algorithm to keep the unconditionally modified energy dissipation law.
* Introducing an adaptive time stepping strategy and restart strategy in order to speed the training process.

The organization of this paper is as follows: In Section 2, we introduce the Evolutionary Deep Operator Neural Network for a given PDE problem. In Section 3, we consider the physics law behind the gradient flow problem and apply the SAV method to obtain the energy dissipation law. We proposed a new architecture for neural network, EDE-DeepONet. In Section 4, we presented two adaptive time stepping strategies, where the second one is called restart in some cases. In Section 5, we generally introduced the architecture of the EDE-DeepONet. In Section 6, we implement our neural network to predict solutions of heat equations, parametric heat equations, and Allen-Cahn equations to show the numerical results.

## 2 Evolutionary Deep Operator Neural Network

Consider a general gradient flow problem,

\[\begin{split}&\frac{\partial\mathbf{u}}{\partial t}+\mathcal{N}_{x}( \mathbf{u})=0\\ &\mathbf{u}(\mathbf{x},0)=\mathbf{f}(\mathbf{x})\end{split} \tag{5}\]

where \(\mathbf{u}\in\mathbf{R}^{l}\), \(\mathcal{N}_{x}(\mathbf{u})\) can be written as a variational derivative of a free energy functional \(E[u(\mathbf{x})]\) bounded from below, \(\mathcal{N}_{x}(\mathbf{u})=\frac{\delta E}{\delta u}.\) The first step is to approximate the initial condition operator with DeepONet.
### Operator learning

For an operator \(\mathcal{G}\), \(\mathcal{G}:\mathbf{u}(\mathbf{x})\mapsto\mathbf{f}(\mathbf{x})\), the data feed into the DeepONet is in the form \((\mathbf{u},y,\mathcal{G}(\mathbf{u})(y))\). It is obtained by the given initial conditions. The branch network takes \([\mathbf{u}(\mathbf{x}_{1}),\mathbf{u}(\mathbf{x}_{2}),\cdots,\mathbf{u}(\mathbf{x}_{m})]^{T}\) as the input, which is the numerical representation of \(\mathbf{u}\), and \([\mathbf{b}_{1},\mathbf{b}_{2},\cdots,\mathbf{b}_{p}]^{T}\in\mathbf{R}^{p\times}\), where \(\mathbf{b}_{k}\in\mathbf{R}^{l}\) for \(k=1,2,\cdots,p\), as outputs. The trunk network takes \(\mathbf{y}\) as the input and \([g_{1},g_{2},\cdots,g_{p}]\in\mathbf{R}^{p}\) as outputs. The Unstacked DeepONet net uses FNN as the baseline model and concatenate the function value at sensor locations and the evaluated point together, i.e. \([\mathbf{u}(\mathbf{x}_{1}),\mathbf{u}(\mathbf{x}_{2}),\cdots,\mathbf{u}(\mathbf{x}_{m}),\mathbf{y}]^{T}\). As the equation in the Universal Approximation Theorem for Operators, we can take the product of \(\mathbf{h}\) and \(t\), then we obtain:

\[\mathcal{G}(\mathbf{u})(\mathbf{x})\approx\sum_{k=1}^{p}g_{k}\mathbf{b}_{k} \tag{6}\]

The activation functions are applied to the trunk net in the last layer. There is no bias in this network. However, according to the theorem 1, the generalization error can be reduced by adding bias. We also give the form with bias \(\mathbf{b}_{0}\):

\[\mathcal{G}(\mathbf{u})(\mathbf{x})\approx\sum_{k=1}^{p}g_{k}\mathbf{b}_{k}+\mathbf{b}_{0} \tag{7}\]

As mentioned before, we assumed the initial condition operator has been trained very well. We are going to find the update rule of the parameters to evolve the neural network.

### The evolution of parameters in the neural network

Denoting the parameters in the branch network as \(W_{1}\) and the parameters in the trunk network as \(W_{2}\), \(W_{1}\) and \(W_{2}\) can be regarded a function of \(t\) since they change in every time step. According to the derivative's chain rule, we have

\[\frac{\partial\mathbf{u}}{\partial t}=\frac{\partial\mathbf{u}}{\partial W_{1}}\frac {\partial W_{1}}{\partial t}+\frac{\partial\mathbf{u}}{\partial W_{2}}\frac{ \partial W_{2}}{\partial t} \tag{8}\]

Since \(\mathbf{u}=\sum_{k=1}^{p}g_{k}\mathbf{b}_{k}=\sum_{k=1}^{p}g_{k}(W_{1}(t))\mathbf{b}_{k}( W_{2}(t))\), then

\[\frac{\partial\mathbf{u}}{\partial t}=\sum_{k=1}^{p}\frac{\partial g_{k}(W_{1}(t) )}{\partial W_{1}}\frac{\partial W_{1}}{\partial t}\mathbf{b}_{k}(W_{2}(t))+\sum_ {k=1}^{p}g_{k}(W_{1}(t))\frac{\partial\mathbf{b}_{k}(W_{2}(t))}{\partial W_{2}} \frac{\partial W_{2}}{\partial t} \tag{9}\]

Our objective is to obtain \(\frac{\partial W_{1}}{\partial t}\) and \(\frac{\partial W_{2}}{\partial t}\), the update rule for parameters. It is equivalent to solve a minimization problem,

\[\left[\frac{\partial W_{1}}{\partial t};\frac{\partial W_{2}}{\partial t} \right]=\text{argmin}\mathcal{J}(\gamma_{1},\gamma_{2}) \tag{10}\]

where

\[\mathcal{J}(\gamma_{1},\gamma_{2})=\frac{1}{2}\left\|\sum_{k=1}^{p}\frac{ \partial g_{k}(W_{1}(t))}{\partial W_{1}}\gamma_{1}\mathbf{b}_{k}(W_{2}(t))+\sum_ {k=1}^{p}g_{k}(W_{1}(t))\frac{\partial\mathbf{b}_{k}(W_{2}(t))}{\partial W_{2}} \gamma_{2}-\mathcal{N}_{\mathbf{x}}(\mathbf{u})\right\|_{2}^{2} \tag{11}\]

In this article, the inner product \((a,b)\) is defined in the integral sense, \((a,b)=\int_{\Omega}a(\mathbf{x})b(\mathbf{x})\,\mathrm{d}\mathbf{x}\) and the \(L_{2}\) norm is defined as \(\left\|a\right\|_{2}^{2}=\int_{\Omega}|a(\mathbf{x})|^{2}\,\mathrm{d}\mathbf{x}\).

The minimization problem can be transformed into a linear system by the first-order optimal condition:

\[\frac{\partial\mathcal{J}}{\partial\gamma_{1}}=\int_{\Omega}\!\! \left(\sum_{k=1}^{p}\frac{\partial g_{k}(W_{1}(t))}{\partial W_{1}}\mathbf{b}_{k}( W_{2}(t))\right)^{T}\left(\gamma_{1}\sum_{k=1}^{p}\frac{\partial g_{k}(W_{1}(t))}{ \partial W_{1}}\mathbf{b}_{k}(W_{2}(t))+\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{ \partial\mathbf{b}_{k}(W_{2}(t))}{\partial W_{2}}\gamma_{2}-\mathcal{N}_{\mathbf{x}}( \mathbf{u})\right)\!\mathrm{d}\mathbf{x}=0 \tag{12}\] \[\frac{\partial\mathcal{J}}{\partial\gamma_{2}}=\int_{\Omega}\!\! \left(\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial\mathbf{b}_{k}(W_{2}(t))}{ \partial W_{2}}\right)^{T}\!\left(\gamma_{1}\sum_{k=1}^{p}\frac{\partial g_{k}(W _{1}(t))}{\partial W_{1}}\mathbf{b}_{k}(W_{2}(t))+\sum_{k=1}^{p}g_{k}(W_{1}(t)) \frac{\partial\mathbf{b}_{k}(W_{2}(t))}{\partial W_{2}}\gamma_{2}-\mathcal{N}_{\bm {x}}(\mathbf{u})\right)\!\mathrm{d}\mathbf{x}=0 \tag{13}\]
In this system, the gradient with respect to \(W_{1}(t)\) and \(W_{2}(t)\) can be computed by automatic differentiation at each time step. By denoting

\[(\mathbf{J_{1}})_{ij_{1}} =\sum_{k=1}^{p}\frac{\partial g_{k}(W_{1}(t))}{\partial W_{1}^{ji} }\mathbf{b}_{k}^{i}(W_{2}(t)) \tag{14}\] \[(\mathbf{J_{2}})_{ij_{2}} =\sum_{k=1}^{p}g_{k}(W_{1}(t))\frac{\partial\mathbf{b}_{k}^{i}(W_{2}( t))}{\partial W_{2}^{ji_{2}}}\] (15) \[(\mathbf{N})_{i} =\mathcal{N}\left(\mathbf{u}_{x}^{i}\right) \tag{16}\]

where \(i=1,2,\cdots,l\), \(j_{1}=1,2,\cdots,N_{\text{para}}^{b}\), \(j_{2}=1,2,\cdots,N_{\text{para}}^{t}\). \(N_{\text{para}}^{b}\) is the number of parameters in Branch net and \(N_{\text{para}}^{t}\) is the number of parameters in Trunk net. \(\mathbf{N}\) is generated by the DeepONet, so it can be evaluated at any spatial point. The above integrals can be approximated by numerical methods:

\[\frac{1}{|\Omega|}\int_{\Omega}\left(\sum_{k=1}^{p}\frac{\partial g _{k}(W_{1}(t))}{\partial W_{1}}\mathbf{b}_{k}(W_{2}(t))\right)^{T}\left(\sum_{k=1} ^{p}\frac{\partial g_{k}(W_{1}(t))}{\partial W_{1}}\mathbf{b}_{k}(W_{2}(t)) \right)\mathrm{d}\mathbf{x} =\lim_{l\to\infty}\frac{1}{l}\mathbf{J_{1}^{T}}\mathbf{J_{1}} \tag{17}\] \[\frac{1}{|\Omega|}\int_{\Omega}\left(\sum_{k=1}^{p}g_{k}(W_{1}(t ))\frac{\partial\mathbf{b}_{k}(W_{2}(t))}{\partial W_{2}}\right)^{T}\left(\sum_{k= 1}^{p}g_{k}(W_{1}(t))\frac{\partial\mathbf{b}_{k}(W_{2}(t))}{\partial W_{2}} \right)\mathrm{d}\mathbf{x} =\lim_{l\to\infty}\frac{1}{l}\mathbf{J_{2}^{T}}\mathbf{J_{2}}\] (18) \[\frac{1}{|\Omega|}\int_{\Omega}\left(\sum_{k=1}^{p}\frac{\partial g _{k}(W_{1}(t))}{\partial W_{1}}\mathbf{b}_{k}(W_{2}(t))\right)^{T}\left(\mathcal{ N}_{\mathbf{x}}(\mathbf{u})\right)\mathrm{d}\mathbf{x} =\lim_{l\to\infty}\frac{1}{l}\mathbf{J_{1}^{T}}\mathbf{N} \tag{19}\]

By denoting \(\gamma_{i}^{opt}\) as optimal values of \(\gamma_{i}\), \(i=1,2\), the objective function can be reduced to

\[\mathbf{J_{1}^{T}}\left(\gamma_{1}^{opt}\mathbf{J_{1}}+\gamma_{2 }^{opt}\mathbf{J_{2}}-\mathbf{N}\right) =0 \tag{20}\] \[\mathbf{J_{2}^{T}}\left(\gamma_{1}^{opt}\mathbf{J_{1}}+\gamma_{2 }^{opt}\mathbf{J_{2}}-\mathbf{N}\right) =0 \tag{21}\]

The feasible solutions of the above equations are the approximated time derivatives of \(W_{1}\) and \(W_{2}\).

\[\frac{dW_{1}}{dt} =\gamma_{1}^{opt} \tag{22}\] \[\frac{dW_{2}}{dt} =\gamma_{2}^{opt} \tag{23}\]

where the initial conditions \(W_{1}^{0}\) and \(W_{2}^{0}\) can be determined by DeepONets for initial condition operators. The two ODEs are the updated rules in the neural networks. The simple way to solve them is the explicit Euler method.

\[\frac{W_{1}^{n+1}-W_{1}^{n}}{\Delta t} =\gamma_{1}^{opt} \tag{24}\] \[\frac{W_{2}^{n+1}-W_{2}^{n}}{\Delta t} =\gamma_{2}^{opt} \tag{25}\]

The neural network can calculate the solution of given PDEs at any time step \(t_{n}\) and spatial point \(\mathbf{x}_{i}\) by weights \(W_{1}^{n}\), \(W_{2}^{n}\), spatial points \(\mathbf{x}\) and initial condition \(\mathbf{u}(\mathbf{x})\).

## 3 Energy Dissipative Evolutionary Deep Operator Neural Network

Let's reconsider the given problem.

\[\frac{\partial\mathbf{u}}{\partial t}+\mathcal{N}_{x}(\mathbf{u})=0 \tag{26}\] \[\mathbf{u}(\mathbf{x},0)=\mathbf{f}(\mathbf{x})\]
where \(\mathbf{u}\in\mathbf{R}^{l}\), \(\mathcal{N}_{\mathbf{x}}(\mathbf{u})\) can be written as a variational derivative of a free energy functional \(E[\mathbf{u}(\mathbf{x})]\) bounded from below, \(\mathcal{N}_{\mathbf{x}}(\mathbf{u})=\frac{\delta E}{\delta\mathbf{u}}\). Taking the inner product with \(\mathcal{N}_{\mathbf{x}}(\mathbf{u})\) of the first equation, we obtain the energy dissipation property

\[\frac{dE[\mathbf{u}(\mathbf{x})]}{dt}=\left(\frac{\delta E}{\delta\mathbf{u}},\frac{ \partial\mathbf{u}}{\partial t}\right)=\left(\mathcal{N}_{\mathbf{x}}(\mathbf{u}),\frac{ \partial\mathbf{u}}{\partial t}\right)=-\left(\mathcal{N}_{\mathbf{x}}(\mathbf{u}), \mathcal{N}_{\mathbf{x}}(\mathbf{u})\right)\leq 0 \tag{27}\]

However, it is usually hard for a numerical algorithm to be efficient as well as energy dissipative. Recently, the SAV approach [33] was introduced to construct numerical schemes which is energy dissipative (with a modified energy), accurate, robust and easy to implement. More precisely, assuming \(E[\mathbf{u}(\mathbf{x})]>0\), it introduces a \(r(t)=\sqrt{E[\mathbf{u}(\mathbf{x},t)]}\), and expands the gradient flow problem as

\[\begin{split}&\frac{\partial\mathbf{u}}{\partial t}=-\frac{r}{\sqrt{E( \mathbf{u})}}\mathcal{N}_{\mathbf{x}}\left(\mathbf{u}\right)\\ & r_{t}=\frac{1}{2\sqrt{E(\mathbf{u})}}\left(\mathcal{N}_{\mathbf{x}} \left(\mathbf{u}\right),\frac{\partial\mathbf{u}}{\partial t}\right)\end{split} \tag{28}\]

With \(r(0)=\sqrt{E[\mathbf{u}(\mathbf{x},t)]}\), the above system has a solution \(r(t)\equiv\sqrt{E[\mathbf{u}(\mathbf{x},t)]}\) and \(\mathbf{u}\) being the solution of the original problem.

### First order scheme

By setting \(\mathbf{u}^{n}=\sum_{k=1}^{p}g_{k}\mathbf{b}_{k}\), a first order scheme can be constructed as

\[\begin{split}&\frac{\mathbf{u}^{n+1}-\mathbf{u}^{n}}{\Delta t}=-\frac{r ^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n})\\ &\frac{r^{n+1}-r^{n}}{\Delta t}=\frac{1}{2\sqrt{E(\mathbf{u}^{n})}} \int_{\Omega}\mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n})\frac{\mathbf{u}^{n+1}-\mathbf{u}^{n}}{ \Delta t}dx.\end{split} \tag{29}\]

This is a coupled system of equations for \((r^{n+1},\mathbf{u}^{n+1})\). But it can be easily decoupled as follows. Plugging the first equation into the second one, we obtain:

\[\frac{r^{n+1}-r^{n}}{\Delta t}=-\frac{r^{n+1}}{2E(\mathbf{u}^{n})}\left\|\mathcal{ N}_{\mathbf{x}}(\mathbf{u}^{n})\right\|^{2}, \tag{30}\]

which implies

\[r^{n+1}=\left(1+\frac{\Delta t}{2E(\mathbf{u}^{n})}\left\|\mathcal{N}_{\mathbf{x}}(\bm {u}^{n})\right\|^{2}\right)^{-1}r^{n} \tag{31}\]

**Theorem 3.1** (Discrete Energy Dissipation Law).: _With the modified energy define above, the scheme is unconditionally energy stable, i.e._

\[(r^{n+1})^{2}-(r^{n})^{2}\leq 0. \tag{32}\]

**Proof 3.1**.: _Taking the inner product of the first equation with \(\frac{r^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n})\) and the second equation with \(2r^{n+1}\)_

\[\begin{split}(r^{n+1})^{2}-(r^{n})^{2}&=2r^{n+1}( r^{n+1}-r^{n})-(r^{n+1}-r^{n})^{2}\\ &=\frac{\Delta r^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\int_{\Omega} \mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n})\frac{\mathbf{u}^{n+1}-\mathbf{u}^{n}}{\Delta t}dx-( r^{n+1}-r^{n})^{2}\\ &=-\left(\frac{r^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\right)^{2}\int_{ \Omega}\mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n})\mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n})dx-(r^{ n+1}-r^{n})^{2}\\ &\leq 0\end{split} \tag{33}\]

In order to maintain the modified energy dissipation law in the evolution neural network, we only need to replace \(\mathcal{N}_{\mathbf{x}}(\mathbf{u})\) by \(\frac{r^{n+1}}{\sqrt{E(\mathbf{u})}}\mathcal{N}_{\mathbf{x}}(\mathbf{u})\) in section 2. The update rule of the neural network is

\[\left[\frac{\partial W_{1}}{\partial t};\,\frac{\partial W_{2}}{\partial t} \right]=\text{argmin}\mathcal{J}(\gamma_{1},\gamma_{2}) \tag{34}\]
where

\[\mathcal{J}(\gamma_{1},\gamma_{2})=\frac{1}{2}\left\|\sum_{k=1}^{p}\frac{\partial g _{k}(W_{1}^{n})}{\partial W_{1}^{n}}\gamma_{1}\mathbf{b}_{k}(W_{2}^{n})+\sum_{k=1}^ {p}g_{k}(W_{1}^{n})\frac{\partial\mathbf{b}_{k}(W_{2}^{n})}{\partial W_{2}^{n}} \gamma_{2}-\frac{r^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\mathcal{N}_{\mathbf{x}}(\mathbf{u}^{n} )\right\|_{2}^{2} \tag{35}\]

The corresponding linear system of the first order optimal condition is

\[\mathbf{J}_{1}^{\rm T}\left(\gamma_{1}^{opt}\mathbf{J}_{1}+\gamma_{2}^{ opt}\mathbf{J}_{2}-\frac{r^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\mathbf{N}\right)=0 \tag{36}\] \[\mathbf{J}_{2}^{\rm T}\left(\gamma_{1}^{opt}\mathbf{J}_{1}+\gamma_{2}^{ opt}\mathbf{J}_{2}-\frac{r^{n+1}}{\sqrt{E(\mathbf{u}^{n})}}\mathbf{N}\right)=0 \tag{37}\]

where

\[(\mathbf{J}_{1})_{ij_{1}} =\sum_{k=1}^{p}\frac{\partial g_{k}(W_{1}^{n})}{\partial W_{1}^{ n,j_{1}}}\mathbf{b}_{k}^{i}(W_{2}^{n}) \tag{38}\] \[(\mathbf{J}_{2})_{ij_{2}} =\sum_{k=1}^{p}g_{k}(W_{1}^{n})\frac{\partial\mathbf{b}_{k}^{i}(W_{2} ^{n})}{\partial W_{2}^{n,j_{2}}}\] (39) \[(\mathbf{N})_{i} =\mathcal{N}\left(\mathbf{u}_{\mathbf{x}}^{i}\right) \tag{40}\]

and \(i=1,2,\cdots,l\), \(j_{1}=1,2,\cdots,N_{\text{para}}^{b}\), \(j_{2}=1,2,\cdots,N_{\text{para}}^{r}\). \(N_{\text{para}}^{b}\) is the number of parameters in Branch net and \(N_{\text{para}}^{t}\) is the number of parameters in Trunk net. After getting \(\gamma_{1}^{opt}\) and \(\gamma_{2}^{opt}\), \(W^{n+1}\) can be obtained by the Forward Euler method as equation (24) and (25).

\[W_{1}^{n+1} =W_{1}^{n}+\gamma_{1}^{opt}\Delta t \tag{41}\] \[W_{2}^{n+1} =W_{2}^{n}+\gamma_{2}^{opt}\Delta t \tag{42}\]

## 4 Adaptive time stepping strategy and Restart strategy

One of the advantages of an unconditionally stable scheme is that the adaptive time step can be utilized. Since the coefficient of \(N_{x}\), \(\frac{r^{n+1}}{\sqrt{E^{n}}}\) should be around 1, by denoting \(\xi^{n+1}=\frac{r^{n+1}}{\sqrt{E^{n}}}\), larger \(\Delta t\) is allowed when \(\xi\) is close to 1 and the smaller \(\Delta t\) is needed when \(\xi\) is far away from 1. Thus, a simple adaptive time-stepping strategy can be described as follows:

```
1. Set the tolerance for \(\xi\) as \(\epsilon_{0}\) and \(\epsilon_{1}\), the initial time step \(\Delta t\), the maximum time step \(\Delta t_{max}\) and the minimum time step \(\Delta t_{min}\)
2. Compute \(u^{n+1}\).
3. Compute \(\xi^{n+1}=\frac{r^{n+1}}{\sqrt{E^{n}}}\).
4. If\(|1-\xi^{n+1}|>\epsilon_{0}\), Then\(\Delta t=\max(\Delta t_{min},\Delta t/2)\); Else if\(|1-\xi^{n+1}|<\epsilon_{1}\), Then\(\Delta t=\min(\Delta t_{max},2\Delta t)\). Go to Step 2.
5. Update time step \(\Delta t\).
```

**Algorithm 1** Adaptive time stepping strategy

Another popular strategy to keep \(r\) approximating the original energy \(E\) is to reset the SAV \(r^{n+1}\) to be \(E^{n+1}\) in some scenarios. The specific algorithm is as following:
```
1. Set the tolerance for \(\epsilon_{0}\), \(\epsilon_{1}\) should be some small tolerance, usually \(10^{-1}\) and \(10^{-3}\). The choices for \(\Delta t_{max}\) and \(\Delta t_{min}\) are quite dependent on \(\Delta t\), usually \(\Delta t_{max}=10^{3}\times\Delta t\) and \(\Delta t_{min}=10^{-3}\times\Delta t\). In Algorithm 2, we usually take \(\epsilon_{2}\) as \(2\times 10^{-2}\).
2. Feed \([u(x_{1}),u(x_{2}),\cdots,u(x_{m})]\) into the branch network and \(y\in Y\) into the trunk network. Denote the output of the DeepONet as \(q\).
3. Update the parameters in the DeepONet by minimizing a cost function, where the cost function can be taken as the mean squared error as \(\frac{1}{|Y|}\sum_{y\in Y}\|\mathcal{G}(u)(y)-q\|^{2}\).
4. Once the DeepONet has been trained well, solve the system of equations of (36) and (37) to obtain \(\left[\frac{\partial W_{1}}{\partial t};\frac{\partial W_{2}}{\partial t}\right]\).
5.The value of \(\left[\frac{\partial W_{1}}{\partial t};\frac{\partial W_{2}}{\partial t}\right]\) can be obtained in the current step. Since the parameters \(W_{1}^{n}\) in the branch network, and \(W_{2}^{n}\) in the trunk network are known, \(W_{1}^{n+1}\) and \(W_{2}^{n+1}\) for the next step can be also obtained by the Forward Euler method or Runge-Kutta method.
6. Repeat step 5 until the final time \(T\), where \(T=t_{0}+s\Delta t\), \(t_{0}\) is the initial time of the given PDE, \(\Delta t\) is the time step in step 5 and \(s\) is the number of repeated times of step 5.
7. Output the solution at time \(T\) in the DeepONet with parameters obtained in step 6.
```

**Algorithm 2** Restart strategy

## 6 Numerical Experiments

In this section, we implement EDE-DeepONet to solve heat equations, parametric heat equations, and Allen-Cahn equations to show its performance and accuracy.

### Example 1: Simple heat equations

To show the accuracy of the EDE-DeepONet, we start with the simple heat equation with different initial conditions since we already have the exact solution. A 1D heat equation system can be described by

\[u_{t}=u_{xx} \tag{43}\] \[u(x,0)=f\] (44) \[u(0,t)=u(2,t)=0 \tag{45}\]

By the method of separation of variables, we can derive the solution to the heat equation. If we set \(f(x)=asin(\pi x)\), the solution is \(u(x,t)=asin(\pi x)e^{-\pi^{2}t}\), where \(a\in[1,2]\). The corresponding energy is \(E(u)=\int_{0}^{2}\frac{1}{2}|u_{x}|^{2}dx\approx\Delta x(\sum_{i=1}^{n}\frac{ 1}{2}|u_{x}(x_{i})|^{2})\). With different parameters \(a\), the above equation describes a kind of PDE. The input data samples can be generated as \((a,x,\mathcal{G}(a)(x))\), where \(\mathcal{G}(a)(x)=a\sin(\pi x)\) for specific \(a\) and \(x\). When generating the initial data samples, we choose 50 points from \([0,2)\) uniformly for x and 50 random values of \(a\) from \([1,2]\). The time step when updating the parameters in the neural network is \(2.5\times 10^{4}\). The number of iteration steps is 400. We compared the different solutions with 4 different \(a\), 1.0, 1.5, 1.8, 2.5 every 100 steps. Although \(a=2.5\) is out of the range of training
data, it still performs well in this model. With the exact solution, we also get the error with different \(a\) as Table 1. The error is defined by \(\frac{1}{N_{x}}\sum_{k=1}^{N_{x}}(u(x_{k})-\hat{u}(x_{k}))^{2}\), where \(N_{x}=51\), \(u\) is the solution obtained by EDE-DeepONet and \(\hat{u}\) is the exact solution. To illustrate the relationship between the modified energy and the original energy, we compare \(r^{2}\) and \(E\) at each step as Figure 2. Both energy are actually disspative in the EDE-DeepONet except when restart strategy applied. The restart strategy is used to keep \(r^{2}\) approaching \(E\). The modified energy is initialized when the restart strategy applied. The restart strategy was triggered on the 370th step since the modified energy and the original energy are offset. After that, they are on the same trajectory again. It is clear that the modified energy approaches the original energy before and after the restart strategy applied. In Figure 3, we give the comparison between the exact solution and the solution obtained by EDE-DeepONet. From this simple heat equation, we show that EDE-DeepONet correctly predicts the solution of the PDE. The most important fact is that EDE-DeepONet can not only predict the solution in the training subset range but also the solution out of the training range. For instance, we take \(a=2.5\) while \(a\in[1,2]\) in the training process. EDE-DeepONet shows good accuracy compared to the exact solution as Figure 3 (a)-(d) and Table 1.
\[u_{t}=cu_{xx} \tag{46}\] \[u(x,0)=sin(\pi x)\] (47) \[u(0,t)=u(2,t)=0 \tag{48}\]

This PDE is more complex than the PDE in Example 1 since the parameter is inside the equation. The traditional numerical scheme needs to be run multiple times to deal with the case with different parameters because they are actually different equations. However, we only need to train the EDE-DeepONet once. The training range of \(c\) is chosen as \([1,2)\). We choose 50 points of \(x\) and \(c\) in the same way as example 1. First, we compared the modified energy with the original energy as Figure 4. The energy is not the same as the first example since the energy depends on the parameter \(c\). We compute the average of the energy with different \(c\) to represent the energy of the system. This case is more complex than the first one, so it needs more restarts during the training. Even though the modified energy oscillates when restart strategy used, it keeps decreasing after each restart. Second, we give the error between the solution obtained by the EDE-DeepONet and the reference solution in Table 2, where the reference solution can be obtained explicitly by variable separation method and the error is defined in the same way as example 1. Third, we give the comparison between our solution and the reference solution in Figure 5. Same as example 1, we give the predicted solution of \(c\notin[1,2]\). All of them show the good accuracy. Hence, EDE-DeepONet can actually solve parametric PDEs.

### Example 3: Allen-Cahn equations

The energy in Examples 1 and 2 is quadratic and the right-hand side of the PDE is linear with respect to \(u\). We are going to show the result for the PDE with more complicated energy. The Allen-Cahn equation is a kind of reaction-diffusion equation. It is derived to describe the process of the phase separation. It was developed to solve a problem in the material science area and has been used to represent the moving interfaces in a phase-field model

Figure 3: The heat equation: The solution with 4 different initial conditions \(f(x)=a\sin(\pi x)\). The curve represents the solution obtained by the EDE-DeepONet, and xxx represents the reference solution. The training parameter \(a\) is in the range of \([1,2)\), so we give three examples in this range. We also present the case out of the range. It also shows accuracy in Figure 3-(d).
in fluid dynamics. The Allen-Cahn equation can be treated as a gradient flow in \(L^{2}\) with some specific energy. We discussed the 1D case and 2D case as follows:

#### 6.3.1 1D case

(a) Various initial conditions:

We start with the simple case, 1D Allen-Cahn equation. It can be described by the following equations:

\[u_{t}=u_{xx}-g(x) \tag{49}\] \[u(x,0)=a\sin\pi x\] (50) \[u(-1,t)=u(1,t)=0 \tag{51}\]

The corresponding Ginzburg-Landau free energy \(E[u]=\int_{0}^{1}\frac{1}{2}|u_{x}|^{2}dx+\int_{x=0}^{x=1}G(u)dx\), where \(G(u)=\frac{1}{4e^{2}}(u^{2}-1)^{2}\) and \(g(u)=G^{\prime}(u)=\frac{1}{e^{2}}u(u^{2}-1)\), \(\epsilon=0.1\). The parameter \(\epsilon\) affects the width of the jump when arriving at the steady state as the Figure 7 (c), (j), (o) and (t). In the EDE-DeepONet, we set \(\Delta t=10^{-4}\), the number of spatial points \(N_{x}\) is 51 and the range of \(a\) is \([0.1,0.5]\). We also compared the modified energy and the original energy as Figure 6. The modified energy can approximate well to the original energy even in a much more complicated form. Then, we compared 4 different solutions with different \(a\in[0.1,0.5]\) obtained by the EDE-DeepONet and the reference solution obtained by the SAV method in traditional numerical computation as Figure 7. The error is shown in Table 3, where error is

\begin{table}
\begin{tabular}{||c|c c c c||} \hline Error & \(T=0.025\) & \(T=0.05\) & \(T=0.075\) & \(T=0.1\) \\ \hline \hline \(c=1.2\) & \(1.30\times 10^{-5}\) & \(1.43\times 10^{-5}\) & \(1.35\times 10^{-5}\) & \(1.20\times 10^{-5}\) \\ \hline \(c=1.5\) & \(1.35\times 10^{-5}\) & \(1.27\times 10^{-5}\) & \(9.80\times 10^{-6}\) & \(7.80\times 10^{-6}\) \\ \hline \(c=1.8\) & \(1.17\times 10^{-5}\) & \(1.03\times 10^{-5}\) & \(7.88\times 10^{-5}\) & \(1.83\times 10^{-5}\) \\ \hline \(c=2.5\) & \(2.20\times 10^{-4}\) & \(1.34\times 10^{-4}\) & \(6.02\times 10^{-5}\) & \(7.08\times 10^{-6}\) \\ \hline \end{tabular}
\end{table}
Table 2: The parametric heat equation: The initial condition of the PDE is \(f(x)=\sin\left(\pi x\right)\). The error is defined by \(\frac{1}{N_{x}}\sum_{k=1}^{N_{x}}(u(x_{0})-\hat{u}(x_{k}))^{2}\), where \(N_{x}=51\), \(u\) is the solution obtained by EDE-DeepONet and \(\hat{u}\) is the exact solution.

Figure 4: The parametric heat equation: The modified energy and original energy when training the network. Each iteration step represents one forward step of the PDE’s numerical solution with \(\Delta t=2.5\times 10^{-4}\). This kind of PDEs is more complicated, so it need more restarts in the training process. The original energy keeps decreasing and the modified energy also shows good approximation of the original energy.
defined in the same way as example 1. \(a=0.6\notin[0.1,0.5)\) shows that EDE-DeepONet can predict the solution well out of the training range. We compared the solution with 4 different initial condition parameter \(a\) every 100 steps until the final time \(T=0.04\) as Figure 7. Each row presents the solution under the same initial condition but with different evolution time \(T\). With this example, it shows that EDE-DeepONet can deal with the PDE with a jump, while it is hard for other neural networks.

(b) Various thickness of the interface:

Heuristically, \(\epsilon\) represents the thickness of the interface in the phase separation process. We are able to obtain a sharp interface when \(\epsilon\to 0\) with evolving in time. Each theoretical and numerical analysis of the limit makes a difference in the purpose of the understanding of the equation, cf. e.g. [34; 35]. We take \(\epsilon\) as a training parameter. The problem can be described as:

\[u_{t}=u_{xx}-\frac{1}{\epsilon^{2}}(u^{3}-u) \tag{52}\] \[u(-1,t)=u(1,t)=0 \tag{53}\]

Figure 5: The parametric heat equation: The solution with 4 different parameters \(c\). The curve represents the solution obtained by the EDE-DeepONet and xxx represents the reference solution. The training parameter \(c\) is in the range of \([1,2)\), so we give 3 examples in this range. We also present the case out of the range in Figure 5-(d).
Since the training sample contains the parameter \(\epsilon\), we can not use the same initial condition as the last example. We use spectral methods for a few steps with initial condition \(u(x,0)=0.4\sin{(\pi x)}\). The training sample is generated based on the numerical solution of \(u_{\epsilon}(x,0.02)\). We randomly select \(50\) different \(\epsilon\) from \([0.1,0.2]\). We set the learning rate as \(\Delta t=10^{-4}\) and apply the adaptive time stepping strategy. We obtain the predicted solution after \(400\) iterations with different \(\epsilon\). The rest setting is the same as the last example. The solution with different \(\epsilon\) is shown in Figure 8. As \(\epsilon\) goes smaller, the interface is sharper. Besides, the range of the training parameter is \((0.1,0.2)\). We are also able to obtain the solution out of the above range. EDE-DeepONet can track the limit of \(\epsilon\) in only one training process, EDE-DeepONet can track the limit of \(\epsilon\) in only one training process, whereas other traditional numerical methods hardly make it.

#### 6.3.2 2D case

The 2D case Allen-Cahn equation is even more complex. The problem can be described as follows:

\[u_{t}=\Delta u-g(u) \tag{54}\] \[u(x,y,0)=asin(\pi x)sin(\pi y)\] (55) \[u(-1,y,t)=u(1,y,t)=u(x,-1,t)=u(x,1,t)=0 \tag{56}\]

The corresponding Ginzburg-Landau free energy \(E[u]=\int_{-1}^{1}\int_{-1}^{1}\frac{1}{2}(|u_{x}|^{2}+|u_{y}|^{2})dxdy+\int_ {-1}^{1}\int_{-1}^{1}G(u)dx\), where \(G(u)=\frac{1}{4e^{2}}(u^{2}-1)^{2}\) and \(g(u)=G^{\prime}(u)=\frac{1}{e}u(u^{2}-1)\). Usually, we take \(\epsilon=0.1\). In the training process, we take \(\Delta t=2\times 10^{-4}\). The number of spatial points is \(51\times 51\) and the number of training parameters \(a\) is \(20\). The way to choose \(a\in(0.1,0.4)\) and \(x\) is the same as in example 1. We first compared the exact solution and the solution obtained by EDE-DeepONet with initial condition \(f(x,y)=0.2sin(\pi x)sin(\pi y)\), where the exact solution is obtained by the traditional SAV method. EDE-DeepONet predicts the solution correctly based on Table 4 and Figure 9. Then in order to show its accuracy, we draw Figure 10 with more parameters. All the examples show good trends to separate. The case \(a=0.4\) is out of the training range, but it still approaches the exact solution.

Figure 6: 1D Allen-Cahn equation: The modified energy and original energy when training the network are shown above. Each iteration step represents one forward step of the PDE’s numerical solution with \(\Delta t=10^{-4}\). The modified energy shows the same trends as the original energy.
## 7 Concluding Remarks

In this paper, we provide a new neural network architecture to solve parametric PDEs with different initial conditions, while maintaining the energy dissipative of dynamic systems. We first introduce the energy dissipative law of dynamic systems to the DeepONet. We also introduce an adaptive time stepping strategy and restart strategy. With our experiments, both above strategies help keep the modified energy approaching the original energy. To avoid much cost of training the DeepONet, we evolve the neural network based on Euler methods. In this article, we adopt the SAV method to solve gradient flow problems. With this successful attempt, more work could be done. For example, we can consider a general Wasserstein gradient flow problem. We are only adopting the basic architecture of the DeepONet. The more advanced architecture is compatible to our work. It may further improve the accuracy of EDE-DeepONet.

## Acknowledgments

SJ and SZ gratefully acknowledge the support of NSF DMS-1720442 and AFOSR FA9550-20-1-0309. GL and ZZ gratefully acknowledge the support of the National Science Foundation (DMS-1555072, DMS-2053746, and DMS

\begin{table}
\begin{tabular}{||c|c c c||} \hline Error & \(T=0.01\) & \(T=0.02\) & \(T=0.03\) \\ \hline \hline \(a=0.15\) & \(1.23\times 10^{-4}\) & \(6.53\times 10^{-4}\) & \(2.75\times 10^{-3}\) \\ \hline \(a=0.2\) & \(2.24\times 10^{-4}\) & \(1.10\times 10^{-3}\) & \(4.04\times 10^{-3}\) \\ \hline \(a=0.3\) & \(4.28\times 10^{-4}\) & \(1.84\times 10^{-3}\) & \(5.76\times 10^{-3}\) \\ \hline \(a=0.35\) & \(5.31\times 10^{-4}\) & \(2.17\times 10^{-3}\) & \(6.25\times 10^{-3}\) \\ \hline \(a=0.4\) & \(6.94\times 10^{-4}\) & \(2.71\times 10^{-3}\) & \(7.22\times 10^{-3}\) \\ \hline \end{tabular}
\end{table}
Table 4: 2D Allen-Cahn equation: The initial condition of the 2D Allen-Cahn equation is \(f(x,y)=a\sin\left(\pi x\right)\sin\left(\pi y\right)\). The error is defined by \(\frac{1}{N_{x}N_{y}}\sum_{k=1}^{N_{x}}\sum_{j=1}^{N_{y}}\sum_{j=1}^{N_{y}}(u(x _{k},y_{j})-\hat{u}(x_{k},y_{j}))^{2}\), where \(N_{x}=N_{y}=51\), \(u\) is the solution obtained by EDE-DeepONet and \(\hat{u}\) is the reference solution.

Figure 7: 1d Allen-Cahn equation: The solution for 1d Allen-Cahn equation with 4 different initial conditions \(f(x)=a\sin\pi x\). The curve represents the solution obtained by our model, and xxx represents the reference solution. We draw the figure for every 100 steps. The range of \(a\) is \([0.1,0.5]\). We also compare the solution with \(a\notin[0.1,0.5]\). All the figures show the trends of the phase separation.
2134209), Brookhaven National Laboratory Subcontract 382247, and U.S. Department of Energy (DOE) Office of Science Advanced Scientific Computing Research program DE-SC0021142 and DE-SC0023161.

Figure 8: 1d Allen-Cahn equation: Solutions with different thickness of the interface at the same final time. The curve represents the solution obtained by EDE-DeepONet. xxx represents the reference solution.
Figure 10: 2D Allen-Cahn equation: The solution of 2D Allen-Cahn equation with 4 different initial conditions \(f(x,y)=a\sin\pi x\sin\pi y\). The training parameter \(a\in[0.1,0.4]\). We draw three figures where \(a\) is in the training range and one figure where \(a\) is out of the training range. All the figures show the phase separation trends according to the reference solution. As \(a\) is further away from the training range, the error tends to be larger.

Figure 9: 2D Allen-Cahn equation: (a)-(d) represents the reference solution of the 2D Allen-Cahn equation with initial condition \(f(x,y)=0.3\sin\left(\pi x\right)\sin\left(\pi y\right)\). (e)-(h) is the solution obtained by the EDE-DeepONet.