【Hackathon 5th No.60】 PhyGeoNet: Physics-Informed Geometry-Adaptive Convolutional Neural Networks for Solving Parameterized Steady-State PDEs on Irregular Domain

docs/zh/examples/phygeonet.md

@@ -5,7 +5,14 @@
     ``` sh
 
     python main.py
+    ```
+=== "模型评估命令"
 
+    ``` sh
+
+    python main.py    mode=eval
+
+    ```
 ## 1. 背景简介
  最近几年，深度学习在很多领域取得了非凡的成就，尤其是计算机视觉和自然语言处理方面，而受启发于深度学习的快速发展，基于深度学习强大的函数逼近能力，神经网络在科学计算领域也取得了成功，现阶段的研究主要分为两大类，一类是将物理信息以及物理限制加入损失函数来对神经网络进行训练, 其代表有 PINN 以及 Deep Retz Net，另一类是通过数据驱动的深度神经网络算子，其代表有 FNO 以及 DeepONet。这些方法都在科学实践中获得了广泛应用，比如天气预测，量子化学，生物工程，以及计算流体等领域。
 ## 2. 问题定义
@@ -18,28 +25,28 @@
 本文使用提出的USCNN模型进行训练。
 ``` py linenums="179"
 --8<--
-examples/hpinns/holography.py:179:188
+examples/phygeonet/heat_equation/main.py:179:188
 --8<--
 ```
 ### 3.2 数据读取
 我们从数据集中读取数据，公开数据集在[AIstudio](https://aistudio.baidu.com/datasetdetail/253292)
 ``` py linenums="166"
 --8<--
-examples/hpinns/holography.py:166:172
+examples/phygeonet/heat_equation/main.py:166:172
 --8<--
 ```
 ### 3.3 输出转化函数构建
 本文为强制边界约束，使用相对应的输出转化函数
 ``` py linenums="215"
 --8<--
-examples/hpinns/holography.py:215:235
+examples/phygeonet/heat_equation/main.py:215:235
 --8<--
 ```
 ### 3.4 约束构建
 构建相对应约束条件，由于边界约束为强制约束，约束条件主要为残差约束
 ``` py linenums="192"
 --8<--
-examples/hpinns/holography.py:192:211
+examples/phygeonet/heat_equation/main.py:192:211
 --8<--
 ```
 
@@ -48,14 +55,14 @@ examples/hpinns/holography.py:192:211
 
 ``` py linenums="254"
 --8<--
-examples/hpinns/holography.py:254:282
+examples/phygeonet/heat_equation/main.py:254:282
 --8<--
 ```
 ### 3.6 优化器构建
 与论文中描述相同，我们使用恒定学习率构造Adam优化器。
 ``` py linenums="190"
 --8<--
-examples/hpinns/holography.py:190:190
+examples/phygeonet/heat_equation/main.py:190:190
 --8<--
 ```
 
@@ -64,15 +71,15 @@ examples/hpinns/holography.py:190:190
 
 ``` py linenums="239"
 --8<--
-examples/epnn/epnn.py:239:247
+examples/phygeonet/heat_equation/main.py:239:247
 --8<--
 ```
 
 最后启动训练即可：
 
 ``` py linenums="249"
 --8<--
-examples/epnn/epnn.py:249:249
+examples/phygeonet/heat_equation/main.pyn.py:249:249
 --8<--
 ```
 
@@ -81,7 +88,6 @@ examples/epnn/epnn.py:249:249
 ``` py linenums="1" title="phygeonet.py"
 --8<--
 examples/phygeonet/heat_equation/main.py
-examples/phygeonet/heat_equation_parameterized_bc/main.py
 --8<--
 ```
 ## 5. 结果展示
@@ -101,6 +107,6 @@ examples/phygeonet/heat_equation_parameterized_bc/main.py
 ## 6. 总结
 本文通过使用调和映射构造坐标变换函数，使得物理信息网络可以在不规则非均匀网格上面进行训练，同时，因为该映射为使用传统方法进行，所以无需训练即可在网络前后嵌入。通过大量实验表明，该网络可以在各种不规则网格问题上表现比SOAT网络突出。
 ## 7. 参考资料
-[NSFnets (Navier-Stokes Flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations](https://arxiv.org/abs/2003.06496)
+[PhyGeoNet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state PDEs on irregular domain](https://www.sciencedirect.com/science/article/pii/S0021999120308536?via%3Dihub)
 
-[Github NSFnets](https://github.com/Alexzihaohu/NSFnets/tree/master)
+[Github PhyGeoNet](https://github.com/Jianxun-Wang/phygeonet/tree/master?tab=readme-ov-file)

examples/phygeonet/heat_equation/conf/conf.yaml

@@ -22,11 +22,10 @@ hydra:
 # general settings
 mode: train # running mode: train/eval
 seed: 66
-date_dir: ./heat_equation.npz
+data_dir: ./data/heat_equation.npz
 output_dir: ${hydra:run.dir}
-boundary_dir: TemplateCase/30/C
 dtheta: 0
-
+epochs: 500
 len_data: 1
 # model settings
 MODEL:
@@ -36,18 +35,16 @@ MODEL:
   h: 0.01
   ny: 19
   nx: 84
-  NvarInput: 2
-  NvarOutput: 1
-  padSingleSide: 1
+  nvar_in: 2
+  nvar_out: 1
+  pad_singleside: 1
 
 # training settings
 TRAIN:
-  epochs: 1000
-  eval_with_no_grad: true
   lr: 1e-3
   batchsize: 1
 
 # evaluation settings
 EVAL:
-  pretrained_model_path: null
+  pretrained_model_path: outputs_phygeonet/latest
   eval_with_no_grad: true

examples/phygeonet/heat_equation/main.py

@@ -1,9 +1,69 @@
+from os import path as osp
+
 import hydra
+import matplotlib
+import matplotlib.pyplot as plt
 import numpy as np
 import paddle
+from matplotlib.collections import PatchCollection
+from matplotlib.patches import Polygon
 from omegaconf import DictConfig
 
 import ppsci
+from ppsci.utils import logger
+
+
+def setAxisLabel(ax, type):
+    if type == "p":
+        ax.set_xlabel(r"$x$")
+        ax.set_ylabel(r"$y$")
+    elif type == "r":
+        ax.set_xlabel(r"$\xi$")
+        ax.set_ylabel(r"$\eta$")
+    else:
+        raise ValueError("The axis type only can be reference or physical")
+
+
+def gen_e2vcg(x):
+    nelemx = x.shape[1] - 1
+    nelemy = x.shape[0] - 1
+    nelem = nelemx * nelemy
+    nnx = x.shape[1]
+    e2vcg = np.zeros([4, nelem])
+    for j in range(nelemy):
+        for i in range(nelemx):
+            e2vcg[:, j * nelemx + i] = np.asarray(
+                [j * nnx + i, j * nnx + i + 1, (j + 1) * nnx + i, (j + 1) * nnx + i + 1]
+            )
+    return e2vcg.astype("int")
+
+
+def visualize2D(ax, x, y, u, colorbarPosition="vertical", colorlimit=None):
+    xdg0 = np.vstack([x.flatten(order="C"), y.flatten(order="C")])
+    udg0 = u.flatten(order="C")
+    idx = np.asarray([0, 1, 3, 2])
+    nelemx = x.shape[1] - 1
+    nelemy = x.shape[0] - 1
+    nelem = nelemx * nelemy
+    e2vcg0 = gen_e2vcg(x)
+    udg_ref = udg0[e2vcg0]
+    cmap = matplotlib.cm.coolwarm
+    polygon_list = []
+    for i in range(nelem):
+        polygon_ = Polygon(xdg0[:, e2vcg0[idx, i]].T)
+        polygon_list.append(polygon_)
+    polygon_ensemble = PatchCollection(polygon_list, cmap=cmap, alpha=1)
+    polygon_ensemble.set_edgecolor("face")
+    polygon_ensemble.set_array(np.mean(udg_ref, axis=0))
+    if colorlimit is None:
+        pass
+    else:
+        polygon_ensemble.set_clim(colorlimit)
+    ax.add_collection(polygon_ensemble)
+    ax.set_xlim(np.min(xdg0[0, :]), np.max(xdg0[0, :]))
+    ax.set_ylim(np.min(xdg0[1, :]), np.max(xdg0[1, :]))
+    cbar = plt.colorbar(polygon_ensemble, orientation=colorbarPosition)
+    return ax, cbar
 
 
 def dfdx(f, dydeta, dydxi, Jinv, h=0.01):
@@ -163,29 +223,17 @@ def train(cfg: DictConfig):
     # initiallizer
     SEED = cfg.seed
     ppsci.utils.misc.set_random_seed(SEED)
-    data = np.load(cfg.date_dir)
+    logger.init_logger("ppsci", osp.join(cfg.output_dir, "train.log"), "info")
+    data = np.load(cfg.data_dir)
     coords = data["coords"]
     jinvs = data["jinvs"]
     dxdxis = data["dxdxis"]
     dydxis = data["dydxis"]
     dxdetas = data["dxdetas"]
     dydetas = data["dydetas"]
-    nx = cfg.MODEL.nx
-    ny = cfg.MODEL.ny
-
-    NvarInput = cfg.MODEL.NvarInput
-    NvarOutput = cfg.MODEL.NvarOutput
-    padSingleSide = cfg.MODEL.padSingleSide
-    model = ppsci.arch.USCNN(
-        cfg.MODEL.input_keys,
-        cfg.MODEL.output_keys,
-        cfg.MODEL.h,
-        nx,
-        ny,
-        NvarInput,
-        NvarOutput,
-        padSingleSide,
-    )
+
+    pad_singleside = cfg.MODEL.pad_singleside
+    model = ppsci.arch.USCNN(**cfg.MODEL)
 
     optimizer = ppsci.optimizer.Adam(cfg.TRAIN.lr)(model)
 
@@ -210,9 +258,9 @@ def train(cfg: DictConfig):
         name="mRes",
     )
 
-    constraint_pde = {sup_constraint_mres.name: sup_constraint_mres}
+    sup_constraint = {sup_constraint_mres.name: sup_constraint_mres}
 
-    def _transform_out(_input, _output, padSingleSide=padSingleSide):
+    def _transform_out(_input, _output, pad_singleside=pad_singleside):
         outputV = _output["outputV"]
         batchSize = outputV.shape[0]
         Jinv = _input["jinvs"]
@@ -221,10 +269,10 @@ def _transform_out(_input, _output, padSingleSide=padSingleSide):
         dxdeta = _input["dxdetas"]
         dydeta = _input["dydetas"]
         for j in range(batchSize):
-            outputV[j, 0, -padSingleSide:, padSingleSide:-padSingleSide] = 0
-            outputV[j, 0, :padSingleSide, padSingleSide:-padSingleSide] = 1
-            outputV[j, 0, padSingleSide:-padSingleSide, -padSingleSide:] = 1
-            outputV[j, 0, padSingleSide:-padSingleSide, 0:padSingleSide] = 1
+            outputV[j, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
+            outputV[j, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
+            outputV[j, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
+            outputV[j, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
             outputV[j, 0, 0, 0] = 0.5 * (outputV[j, 0, 0, 1] + outputV[j, 0, 1, 0])
             outputV[j, 0, 0, -1] = 0.5 * (outputV[j, 0, 0, -2] + outputV[j, 0, 1, -1])
         dvdx = dfdx(outputV, dydeta, dydxi, Jinv)
@@ -238,12 +286,11 @@ def _transform_out(_input, _output, padSingleSide=padSingleSide):
     output_dir = cfg.output_dir
     solver = ppsci.solver.Solver(
         model,
-        constraint_pde,
+        sup_constraint,
         output_dir,
         optimizer,
-        epochs=cfg.TRAIN.epochs,
+        epochs=cfg.epochs,
         iters_per_epoch=coords.shape[0],
-        eval_with_no_grad=cfg.TRAIN.eval_with_no_grad,
     )
 
     solver.train()
@@ -255,13 +302,13 @@ def evaluate(cfg: DictConfig):
     SEED = cfg.seed
     ppsci.utils.misc.set_random_seed(SEED)
 
-    data = np.load(cfg.date_dir)
+    data = np.load(cfg.data_dir)
     coords = data["coords"]
 
     OFV_sb = paddle.to_tensor(data["OFV_sb"])
 
     ## create model
-    padSingleSide = 1
+    pad_singleside = 1
     model = ppsci.arch.USCNN(**cfg.MODEL)
     model.register_output_transform(None)
     solver = ppsci.solver.Solver(
@@ -270,17 +317,47 @@ def evaluate(cfg: DictConfig):
     )
     outputV = solver.predict({"coords": paddle.to_tensor(coords)})
     outputV = outputV["outputV"]
-    outputV[0, 0, -padSingleSide:, padSingleSide:-padSingleSide] = 0
-    outputV[0, 0, :padSingleSide, padSingleSide:-padSingleSide] = 1
-    outputV[0, 0, padSingleSide:-padSingleSide, -padSingleSide:] = 1
-    outputV[0, 0, padSingleSide:-padSingleSide, 0:padSingleSide] = 1
+    outputV[0, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
+    outputV[0, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
+    outputV[0, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
+    outputV[0, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
     outputV[0, 0, 0, 0] = 0.5 * (outputV[0, 0, 0, 1] + outputV[0, 0, 1, 0])
     outputV[0, 0, 0, -1] = 0.5 * (outputV[0, 0, 0, -2] + outputV[0, 0, 1, -1])
     CNNVNumpy = outputV[0, 0, :, :]
     ev = paddle.sqrt(
         paddle.mean((OFV_sb - CNNVNumpy) ** 2) / paddle.mean(OFV_sb**2)
     ).item()
     print(ev)
+    outputV = outputV.numpy()
+    OFV_sb = OFV_sb.numpy()
+    fig1 = plt.figure()
+    ax = plt.subplot(1, 2, 1)
+    visualize2D(
+        ax,
+        coords[0, 0, 1:-1, 1:-1],
+        coords[0, 1, 1:-1, 1:-1],
+        outputV[0, 0, 1:-1, 1:-1],
+        "horizontal",
+        [0, 1],
+    )
+    setAxisLabel(ax, "p")
+    ax.set_title("CNN " + r"$T$")
+    ax.set_aspect("equal")
+    ax = plt.subplot(1, 2, 2)
+    visualize2D(
+        ax,
+        coords[0, 0, 1:-1, 1:-1],
+        coords[0, 1, 1:-1, 1:-1],
+        OFV_sb[1:-1, 1:-1],
+        "horizontal",
+        [0, 1],
+    )
+    setAxisLabel(ax, "p")
+    ax.set_aspect("equal")
+    ax.set_title("FV " + r"$T$")
+    fig1.tight_layout(pad=1)
+    fig1.savefig("T.pdf", bbox_inches="tight")
+    plt.close(fig1)
 
 
 @hydra.main(version_base=None, config_path="./conf", config_name="conf.yaml")

examples/phygeonet/heat_equation_parameterized_BC/conf/conf.yaml

@@ -20,15 +20,16 @@ hydra:
     subdir: ./
 
 # general settings
-mode: train # running mode: train/eval
+mode: eval # running mode: train/eval
 seed: 66
-date_dir: heat_equation_bc.npz
+data_dir: ./data/heat_equation_bc.npz
+test_data_dir: ./data/heat_equation_bc_test.npz
 output_dir: ${hydra:run.dir}
 r: 0.5
 R: 1
 dtheta: 0
 len_data: 2
-casename: ["TemplateCase0", "TemplateCase6"]
+epochs: 2000
 # model settings
 MODEL:
   input_keys: [ 'coords' ]
@@ -37,18 +38,16 @@ MODEL:
   h: 0.01
   ny: 276
   nx: 49
-  nVarIn: 1
-  nVarOut: 1
-  padSingleSide: 1
+  nvar_in: 1
+  nvar_out: 1
+  pad_singleside: 1
 
 # training settings
 TRAIN:
-  epochs: 2000
-  eval_with_no_grad: true
   lr: 1e-3
   batchsize: 2
 
 # evaluation settings
 EVAL:
-  pretrained_model_path: ./checkpoints/latest
+  pretrained_model_path: outputs_phygeonet/2023-12-19/22-59-16/checkpoints/latest
   eval_with_no_grad: true