Add documentation for diffusion-reaction demo (lululxvi#1012)

docs/demos/pinn_forward.rst

@@ -44,12 +44,12 @@ Time-dependent PDEs
    pinn_forward/diffusion.1d
    pinn_forward/diffusion.1d.exactBC
    pinn_forward/diffusion.1d.resample
+   pinn_forward/diffusion.reaction
    pinn_forward/burgers.rar
    pinn_forward/allen.cahn
    pinn_forward/klein.gordon
 
 - `Heat equation with training points resampling <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/heat_resample.py>`_
-- `Diffusion-reaction equation <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/diffusion_reaction.py>`_
 - `Beltrami flow <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/Beltrami_flow.py>`_
 - `Kovasznay flow <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/Kovasznay_flow.py>`_
 - `Wave propagation with spatio-temporal multi-scale Fourier feature architecture <https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/wave_1d.py>`_

docs/demos/pinn_forward/diffusion.reaction.rst

@@ -0,0 +1,180 @@
+Diffusion-reaction equation
+===========================
+
+Problem setup
+--------------
+
+We will solve the following 1D diffusion-reaction equation:
+
+.. math:: \frac{\partial y}{\partial t} = d \frac{\partial^2 y}{\partial x^2} + e^{-t} (3\frac{\sin{2x}}{2} + \frac{8\sin{3x}}{3} + \frac{15\sin{4x}}{4} + \frac{63\sin{8x}}{8})
+
+with the initial condition
+
+.. math:: y(x, 0) = \sin{x} + \frac{\sin{2x}}{2} + \frac{\sin{3x}}{3} + \frac{\sin{4x}}{4} + \frac{\sin{8x}}{8}, \quad x \in [-\pi, \pi]
+
+and the Dirichlet boundary condition
+
+.. math:: y(t, -\pi) = y(t, \pi) = 0, \quad t \in [0, 1]
+
+We also specify the following parameters for the equation:
+
+.. math:: d = 1 
+
+The exact solution is
+
+.. math:: y(x, t) = e^{-t}( \sin{x} + \frac{\sin{2x}}{2} + \frac{\sin{3x}}{3} + \frac{\sin{4x}}{4} + \frac{\sin{8x}}{8})
+
+Implementation
+--------------
+
+This description goes through the implementation of a solver for the above described diffusion-reaction equation step-by-step.
+
+First, DeepXDE, Numpy and Tensorflow libraries are imported:
+
+.. code-block:: python
+
+    import DeepXDE as dde
+    import numpy as np
+    import tensorflow as tf
+
+Next, we express the PDE residual of the diffusion-reaction equation:
+
+.. code-block:: python
+
+    def pde(x, y):
+        dy_t = dde.grad.jacobian(y, x, i=0, j=1)
+        dy_xx = dde.grad.hessian(y, x, i=0, j=0)
+        d = 1
+        # Backend tensorflow.compact.v1 or tensorflow
+        return (
+            dy_t
+            - d * dy_xx
+            - tf.exp(-x[:, 1:])
+            * (
+                3 * tf.sin(2 * x[:, 0:1]) / 2
+                + 8 * tf.sin(3 * x[:, 0:1]) / 3
+                + 15 * tf.sin(4 * x[:, 0:1]) / 4
+                + 63 * tf.sin(8 * x[:, 0:1]) / 8
+            )
+        )
+
+If the backend is Pytorch, the residual should look like this:
+
+.. code-block:: python
+
+    def pde(x, y):
+        dy_t = dde.grad.jacobian(y, x, i=0, j=1)
+        dy_xx = dde.grad.hessian(y, x, i=0, j=0)
+        d = 1
+        # Backend pytorch
+        return (
+            dy_t
+            - d * dy_xx
+            - torch.exp(-x[:, 1:])
+            * (
+                3 * torch.sin(2 * x[:, 0:1]) / 2
+                + 8 * torch.sin(3 * x[:, 0:1]) / 3
+                + 15 * torch.sin(4 * x[:, 0:1]) / 4
+                + 63 * torch.sin(8 * x[:, 0:1]) / 8
+            )
+        )
+
+The first argument to ``pde`` is the 2 dimensional vector where the first component(``x[:, 0]``) is the :math:`x`-coordinate, and the second component(``x[:, 1]``) is the :math:`t`-coordinate. The second argument is the network output, i.e., the solution :math:`u(x)`, but here we use ``y`` as the name of the variable.
+
+Then we define the solution to the PDE:
+
+.. code-block:: python
+
+    def func(x):
+        return np.exp(-x[:, 1:]) * (
+            np.sin(x[:, 0:1])
+            + np.sin(2 * x[:, 0:1]) / 2
+            + np.sin(3 * x[:, 0:1]) / 3
+            + np.sin(4 * x[:, 0:1]) / 4
+            + np.sin(8 * x[:, 0:1]) / 8
+        )
+
+
+Now we can define a computational geometry and time domain. We can use a built-in class ``Interval`` and ``TimeDomain`` and we combine both the domains using ``GeometryXTime`` as follows
+
+.. code-block:: python
+
+    geom = dde.geometry.Interval(-1, 1)
+    timedomain = dde.geometry.TimeDomain(0, 0.99)
+    geomtime = dde.geometry.GeometryXTime(geom, timedomain)
+
+Now, we have specified the geometry and the PDE residual. We then define the ``TimePDE`` problem as
+
+.. code-block:: python
+
+    data = dde.data.TimePDE(
+        geomtime, pde, [], num_domain=320, solution=func, num_test=80000
+    )
+
+The number 320 is the number of training residual points sampled inside the domain, and the number 80000 is the number of points sampled inside
+the domain for testing the PDE loss.
+
+Next, we choose the network. Here, we use a fully connected neural network of depth 7 (i.e., 6 hidden layers) and width 30:
+
+.. code-block:: python
+
+    layer_size = [2] + [30] * 6 + [1]
+    activation = "tanh"
+    initializer = "Glorot uniform"
+    net = dde.nn.FNN(layer_size, activation, initializer)
+
+
+Then we construct a function that satisfies both the initial and the boundary conditions to tansform the network output.
+In this case, :math:`t(\pi^2 - x^2)y + \sin{x} + \frac{\sin{2x}}{2} + \frac{\sin{3x}}{3} + \frac{\sin{4x}}{4} + \frac{\sin{8x}}{8}` is used.
+If :math:`t` is equal to 0, the initial condition is recovered. When :math:`x` is equal to :math:`-\pi` or :math:`\pi`, the boundary condition is recovered.
+Hence the initial and boundary conditions are both hard conditions.
+
+.. code-block:: python
+
+    def output_transform(x, y):
+        return (
+            x[:, 1:2] * (np.pi ** 2 - x[:, 0:1] ** 2) * y
+            + tf.sin(x[:, 0:1])
+            + tf.sin(2 * x[:, 0:1]) / 2
+            + tf.sin(3 * x[:, 0:1]) / 3
+            + tf.sin(4 * x[:, 0:1]) / 4
+            + tf.sin(8 * x[:, 0:1]) / 8
+        )
+    
+    net.apply_output_transform(output_transform)
+
+If the backend is Pytorch, the code should look like this:
+
+.. code-block:: python
+
+    def output_transform(x, y):
+        return (
+            x[:, 1:2] * (np.pi ** 2 - x[:, 0:1] ** 2) * y
+            + torch.sin(x[:, 0:1])
+            + torch.sin(2 * x[:, 0:1]) / 2
+            + torch.sin(3 * x[:, 0:1]) / 3
+            + torch.sin(4 * x[:, 0:1]) / 4
+            + torch.sin(8 * x[:, 0:1]) / 8
+        )
+    
+    net.apply_output_transform(output_transform)
+
+Now, we have the PDE problem and the network. We build a ``Model`` and choose the optimizer and learning rate. We then train the model for 20000 iterations.
+
+.. code-block:: python
+
+    model = dde.Model(data, net)
+    model.compile("adam", lr=1e-3, metrics=["l2 relative error"])
+    losshistory, train_state = model.train(iterations=20000)
+
+We also save and plot the best trained result and loss history.
+
+.. code-block:: python
+
+    dde.saveplot(losshistory, train_state, issave=True, isplot=True)
+
+Complete code
+-------------
+
+.. literalinclude:: ../../../examples/pinn_forward/diffusion_reaction.py
+  :language: python