赛题十五：PIRBN 子目录 1D_sine_function

jointContribution/PIRBN/README.md

@@ -0,0 +1,34 @@
+# Physics-informed radial basis network (PIRBN)
+
+This repository provides numerical examples of the **physics-informed radial basis network** (**PIRBN**).
+
+Physics-informed neural network (PINN) has recently gained increasing interest in computational  mechanics.
+
+This work starts from studying the training dynamics of PINNs via the nerual tangent kernel (NTK) theory. Based on numerical experiments, we found:
+
+- PINNs tend to be a **local approximator** during the training
+- For PINNs who fail to be a local apprixmator, the physics-informed loss can be hardly minimised through training
+
+Inspired by findings, we proposed the PIRBN, which can exhibit the local property intuitively. It has been demonstrated that the NTK theory is applicable for PIRBN. Besides, other PINN techniques can be directly migrated to PIRBNs.
+
+Numerical examples include:
+
+ - 1D sine funtion (**Eq. 1** in the manuscript)
+
+      **PDE**: $\frac{\partial^2 }{\partial x^2}u(x)-4\mu^2\pi^2 sin(2\mu\pi(x))=0, x\in[0,1]$
+
+      **BC**:  $u(0)=u(1)=0.$
+
+For more details in terms of mathematical proofs and numerical examples, please refer to our paper.
+
+# Link
+
+<https://doi.org/10.1016/j.cma.2023.116290>
+
+<https://github.com/JinshuaiBai/PIRBN>
+
+# Enviornmental settings
+
+```
+pip install -r requirements.txt
+```

jointContribution/PIRBN/analytical_solution.py

@@ -0,0 +1,72 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import paddle
+import scipy.io
+
+
+def output_fig(train_obj, mu, b):
+    plt.figure(figsize=(15, 9))
+    rbn = train_obj.pirbn.rbn
+
+    target_dir = os.path.join(os.path.dirname(__file__), "../target")
+    if not os.path.exists(target_dir):
+        os.mkdir(target_dir)
+
+    # Comparisons between the PINN predictions and the ground truth.
+    plt.subplot(2, 3, 1)
+    ns = 1001
+    dx = 1 / (ns - 1)
+    xy = np.zeros((ns, 1)).astype(np.float32)
+    for i in range(0, ns):
+        xy[i, 0] = i * dx
+    y = rbn(paddle.to_tensor(xy))
+    y = y.numpy()
+    y_true = np.sin(2 * mu * np.pi * xy)
+    plt.plot(xy, y_true)
+    plt.plot(xy, y, linestyle="--")
+    plt.legend(["ground truth", "PINN"])
+
+    # Point-wise absolute error plot.
+    plt.subplot(2, 3, 2)
+    plt.plot(xy, np.abs(y_true - y))
+    plt.ylim(top=1e-3)
+    plt.legend(["Absolute Error"])
+
+    # Loss history of the PINN during the training process.
+    plt.subplot(2, 3, 3)
+    his_l1 = train_obj.his_l1
+    x = range(len(his_l1))
+    plt.yscale("log")
+    plt.plot(x, his_l1)
+    plt.plot(x, train_obj.his_l2)
+    plt.legend(["Lg", "Lb"])
+
+    # Visualise NTK after initialisation, The normalised Kg at 0th iteration.
+    plt.subplot(2, 3, 4)
+    jac = train_obj.ntk_list[0]
+    a = np.dot(jac, np.transpose(jac))
+    plt.imshow(a / (np.max(abs(a))), cmap="bwr", vmax=1, vmin=-1)
+    plt.colorbar()
+
+    # Visualise NTK after training, The normalised Kg at 2000th iteration.
+    plt.subplot(2, 3, 5)
+    if 2000 in train_obj.ntk_list:
+        jac = train_obj.ntk_list[2000]
+        a = np.dot(jac, np.transpose(jac))
+        plt.imshow(a / (np.max(abs(a))), cmap="bwr", vmax=1, vmin=-1)
+        plt.colorbar()
+
+    # The normalised Kg at 20000th iteration.
+    plt.subplot(2, 3, 6)
+    if 20000 in train_obj.ntk_list:
+        jac = train_obj.ntk_list[20000]
+        a = np.dot(jac, np.transpose(jac))
+        plt.imshow(a / (np.max(abs(a))), cmap="bwr", vmax=1, vmin=-1)
+        plt.colorbar()
+
+    plt.savefig(os.path.join(target_dir, f"sine_function_{mu}_{b}.png"))
+
+    # Save data
+    scipy.io.savemat(os.path.join(target_dir, "out.mat"), {"NTK": a, "x": xy, "y": y})

jointContribution/PIRBN/main.py

@@ -0,0 +1,38 @@
+import analytical_solution
+import numpy as np
+import pirbn
+import rbn_net
+import train
+
+# Define mu
+mu = 4
+
+# Define the number of sample points
+ns = 51
+
+# Define the sample points' interval
+dx = 1.0 / (ns - 1)
+
+# Initialise sample points' coordinates
+xy = np.zeros((ns, 1)).astype(np.float32)
+for i in range(0, ns):
+    xy[i, 0] = i * dx
+xy_b = np.array([[0.0], [1.0]])
+
+x = [xy, xy_b]
+y = [-4 * mu**2 * np.pi**2 * np.sin(2 * mu * np.pi * xy)]
+
+# Set up radial basis network
+n_in = 1
+n_out = 1
+n_neu = 61
+b = 10.0
+c = [-0.1, 1.1]
+
+# Set up PIRBN
+rbn = rbn_net.RBN_Net(n_in, n_out, n_neu, b, c)
+train_obj = train.Trainer(pirbn.PIRBN(rbn), x, y, learning_rate=0.001, maxiter=20001)
+train_obj.fit()
+
+# Visualise results
+analytical_solution.output_fig(train_obj, mu, b)

jointContribution/PIRBN/pirbn.py

@@ -0,0 +1,85 @@
+import paddle
+
+
+class Dif(paddle.nn.Layer):
+    """This function is to initialise for differential operator.
+
+    Args:
+        rbn (model): The Feedforward Neural Network.
+    """
+
+    def __init__(self, rbn, **kwargs):
+        super().__init__(**kwargs)
+        self.rbn = rbn
+
+    def forward(self, x):
+        """This function is to calculate the differential terms.
+
+        Args:
+            x (Tensor): The coordinate array
+
+        Returns:
+            Tuple[Tensor, Tensor]: The first-order derivative of the u with respect to the x; The second-order derivative of the u with respect to the x.
+        """
+        x.stop_gradient = False
+        # Obtain the output from the RBN
+        u = self.rbn(x)
+        # Obtain the first-order derivative of the output with respect to the input
+        u_x = paddle.grad(u, x, retain_graph=True, create_graph=True)[0]
+        # Obtain the second-order derivative of the output with respect to the input
+        u_xx = paddle.grad(u_x, x, retain_graph=True, create_graph=True)[0]
+        return u_x, u_xx
+
+
+class PIRBN(paddle.nn.Layer):
+    def __init__(self, rbn):
+        super().__init__()
+        self.rbn = rbn
+
+    def forward(self, input_data):
+        xy, xy_b = input_data
+        # initialize the differential operators
+        Dif_u = Dif(self.rbn)
+        u_b = self.rbn(xy_b)
+
+        # obtain partial derivatives of u with respect to x
+        _, u_xx = Dif_u(xy)
+
+        return [u_xx, u_b]
+
+    def cal_ntk(self, x):
+        # Formula (4), Page5, \gamma variable
+        gamma_g = 0.0
+        gamma_b = 0.0
+        n_neu = self.rbn.n_neu
+
+        # in-domain
+        n1 = x[0].shape[0]
+        for i in range(n1):
+            temp_x = [x[0][i, ...].unsqueeze(0), paddle.to_tensor([[0.0]])]
+            y = self.forward(temp_x)
+            l1t = paddle.grad(y[0], self.parameters())
+            for j in l1t:
+                gamma_g = gamma_g + paddle.sum(j**2) / n1
+            temp = paddle.concat((l1t[0], l1t[1].reshape((1, n_neu))), axis=1)
+            if i == 0:
+                # Fig.1, Page8, Kg variable
+                Kg = temp
+            else:
+                Kg = paddle.concat((Kg, temp), axis=0)
+
+        # bound
+        n2 = x[1].shape[0]
+        for i in range(n2):
+            temp_x = [paddle.to_tensor([[0.0]]), x[1][i, ...].unsqueeze(0)]
+            y = self.forward(temp_x)
+            l1t = paddle.grad(y[1], self.parameters())
+            for j in l1t:
+                gamma_b = gamma_b + paddle.sum(j**2) / n2
+
+        # calculate adapt factors
+        temp = gamma_g + gamma_b
+        gamma_g = temp / gamma_g
+        gamma_b = temp / gamma_b
+
+        return gamma_g, gamma_b, Kg

jointContribution/PIRBN/rbn_net.py

@@ -0,0 +1,88 @@
+import math
+
+import numpy as np
+import paddle
+
+
+class RBN_Net(paddle.nn.Layer):
+    """This class is to build a radial basis network (RBN).
+
+    Args:
+        n_in (int): Number of input of the RBN.
+        n_out (int): Number of output of the RBN.
+        n_neu (int): Number of neurons in the hidden layer.
+        b (List[float32]|float32): Initial value for hyperparameter b.
+        c (List[float32]): Initial value for hyperparameter c.
+    """
+
+    def __init__(self, n_in, n_out, n_neu, b, c):
+        super().__init__()
+        self.n_in = n_in
+        self.n_out = n_out
+        self.n_neu = n_neu
+        self.b = paddle.to_tensor(b)
+        self.c = paddle.to_tensor(c)
+
+        self.layer1 = RBF_layer1(self.n_neu, self.c, n_in)
+        # LeCun normal
+        std = math.sqrt(1 / self.n_neu)
+        self.linear = paddle.nn.Linear(
+            self.n_neu,
+            self.n_out,
+            weight_attr=paddle.ParamAttr(
+                initializer=paddle.nn.initializer.Normal(mean=0.0, std=std)
+            ),
+            bias_attr=False,
+        )
+        self.ini_ab()
+
+    def forward(self, x):
+        temp = self.layer1(x)
+        y = self.linear(temp)
+        return y
+
+    def ini_ab(self):
+        b = np.ones((1, self.n_neu)) * self.b
+        self.layer1.b = self.layer1.create_parameter(
+            (1, self.n_neu), default_initializer=paddle.nn.initializer.Assign(b)
+        )
+
+
+class RBF_layer1(paddle.nn.Layer):
+    """This class is to create the hidden layer of a radial basis network.
+
+    Args:
+        n_neu (int): Number of neurons in the hidden layer.
+        c (List[float32]): Initial value for hyperparameter b.
+        input_shape_last (int): Last item of input shape.
+    """
+
+    def __init__(self, n_neu, c, input_shape_last):
+        super(RBF_layer1, self).__init__()
+        self.n_neu = n_neu
+        self.c = c
+        self.b = self.create_parameter(
+            shape=(input_shape_last, self.n_neu),
+            dtype=paddle.get_default_dtype(),
+            # Convert from tensorflow tf.random_normal_initializer(), default value mean=0.0, std=0.05
+            default_initializer=paddle.nn.initializer.Normal(mean=0.0, std=0.05),
+        )
+
+    def forward(self, inputs):  # Defines the computation from inputs to outputs
+        temp_x = paddle.matmul(inputs, paddle.ones((1, self.n_neu)))
+        x0 = (
+            paddle.reshape(
+                paddle.arange(self.n_neu, dtype=paddle.get_default_dtype()),
+                (1, self.n_neu),
+            )
+            * (self.c[1] - self.c[0])
+            / (self.n_neu - 1)
+            + self.c[0]
+        )
+        x_new = temp_x - x0
+        return self.rbf_activate(x_new)
+
+    # activation function
+    def rbf_activate(self, input):
+        s = self.b * self.b
+        return paddle.exp(-(input * input) * s)

jointContribution/PIRBN/train.py

@@ -0,0 +1,57 @@
+import paddle
+
+# Used to calculate the second-order derivatives
+paddle.framework.core.set_prim_eager_enabled(True)
+
+
+class Trainer:
+    def __init__(self, pirbn, x_train, y_train, learning_rate=0.001, maxiter=10000):
+        # set attributes
+        self.pirbn = pirbn
+
+        self.learning_rate = learning_rate
+        self.x_train = [
+            paddle.to_tensor(x, dtype=paddle.get_default_dtype()) for x in x_train
+        ]
+        self.y_train = paddle.to_tensor(y_train, dtype=paddle.get_default_dtype())
+        self.maxiter = maxiter
+        self.his_l1 = []
+        self.his_l2 = []
+        self.iter = 0
+        self.a_g = paddle.to_tensor(1.0)
+        self.a_b = paddle.to_tensor(1.0)
+        self.optimizer = paddle.optimizer.Adam(
+            learning_rate=0.001, parameters=self.pirbn.parameters()
+        )
+        self.ntk_list = {}
+
+    def Loss(self, x, y, a_g, a_b):
+        tmp = self.pirbn(x)
+        l1 = 0.5 * paddle.mean(paddle.square(tmp[0] - y[0]))
+        l2 = 0.5 * paddle.mean(paddle.square(tmp[1]))
+        loss = l1 * a_g + l2 * a_b
+        return loss, l1, l2
+
+    def evaluate(self):
+        # compute loss
+        loss, l1, l2 = self.Loss(self.x_train, self.y_train, self.a_g, self.a_b)
+        l1_numpy = float(l1)
+        l2_numpy = float(l2)
+        self.his_l1.append(l1_numpy)
+        self.his_l2.append(l2_numpy)
+        if self.iter % 200 == 0:
+            self.a_g, self.a_b, _ = self.pirbn.cal_ntk(self.x_train)
+            print("\ta_g =", float(self.a_g), "\ta_b =", float(self.a_b))
+            print("Iter: ", self.iter, "\tL1 =", l1_numpy, "\tL2 =", l2_numpy)
+        self.iter = self.iter + 1
+        return loss
+
+    def fit(self):
+        for i in range(0, self.maxiter):
+            if i in [0, 2000, 20000]:
+                self.ntk_list[i] = self.pirbn.cal_ntk(self.x_train)[2].numpy()
+            loss = self.evaluate()
+            loss.backward()
+            self.optimizer.step()
+            self.optimizer.clear_grad()
+        return loss, [self.his_l1, self.his_l2]