Skip to content

PhyGeoNet

AI Studio Quick Experience

# heat_equation
# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation.npz -P ./data/

# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation.npz --create-dirs -o ./data/heat_equation.npz

python heat_equation.py

# heat_equation_bc
# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc.npz -P ./data/
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc_test.npz -P ./data/

# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc.npz --create-dirs -o ./data/heat_equation.npz
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc_test.npz --create-dirs -o ./data/heat_equation.npz

python heat_equation_with_bc.py

# heat_equation
# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation.npz -P ./data/

# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation.npz --create-dirs -o ./data/heat_equation.npz

python heat_equation.py mode=eval EVAL.pretrained_model_path=https://paddle-org.bj.bcebos.com/paddlescience/models/PhyGeoNet/heat_equation_pretrain.pdparams

# heat_equation_bc
# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc.npz -P ./data/
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc_test.npz -P ./data/

# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc.npz --create-dirs -o ./data/heat_equation.npz
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc_test.npz --create-dirs -o ./data/heat_equation.npz

python heat_equation_with_bc.py mode=eval EVAL.pretrained_model_path=https://paddle-org.bj.bcebos.com/paddlescience/models/PhyGeoNet/heat_equation_bc_pretrain.pdparams

# heat_equation
python heat_equation.py mode=export

# heat_equation_bc
python heat_equation_with_bc.py mode=export

# heat_equation
# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation.npz -P ./data/

# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation.npz --create-dirs -o ./data/heat_equation.npz

python heat_equation.py mode=infer

# heat_equation_bc
# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc.npz -P ./data/
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc_test.npz -P ./data/

# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc.npz --create-dirs -o ./data/heat_equation.npz
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyGeoNet/heat_equation_bc_test.npz --create-dirs -o ./data/heat_equation.npz

python heat_equation_with_bc.py mode=infer
Model mRes ev
heat_equation_pretrain.pdparams 0.815 0.095
heat_equation_bc_pretrain.pdparams 992 0.31

1. Background Introduction

In recent years, deep learning has achieved remarkable achievements in many fields, especially in computer vision and natural language processing. Inspired by the rapid development of deep learning and based on the powerful function approximation ability of deep learning, neural networks have also achieved success in the field of scientific computing. Current research is mainly divided into two categories. One is to add physical information and physical constraints to the loss function to train neural networks, represented by PINN and Deep Ritz Net. The other is data-driven deep neural network operators, represented by FNO and DeepONet. These methods have been widely used in scientific practice, such as weather forecasting, quantum chemistry, biological engineering, and computational fluid dynamics. Due to the parameter sharing nature of convolutional neural networks, they can learn large-scale spatiotemporal domains, so they have received more and more attention.

2. Problem Definition

In actual scientific computing problems, the solution domain of many partial differential equations has complex boundaries and is non-uniform. Existing neural networks often target solution domains with regular boundaries and uniform grids, so they have no practical application effect.

Aiming at the problem that physical information neural networks perform poorly on complex boundary non-uniform grid solution domains, this paper proposes a method to transform irregular boundary non-uniform grids into regular boundary uniform grids through coordinate transformation. In addition, this paper uses the above advantages of convolutional neural networks after becoming uniform grids, and proposes corresponding physical information convolutional neural networks.

3. Problem Solving

To save space, heat equation will be used as an example to explain how to implement it using PaddleScience.

3.1 Model Construction

This case uses the proposed USCNN model for training. The construction method of this model is shown below.

23

Among them, the parameters required to build the model can be obtained from the corresponding configuration file.

34
35
36
37
38
39
40
41
42
43
# model settings
MODEL:
  input_keys: ["coords"]
  output_keys: ["output_v"]
  hidden_size: [16, 32, 16]
  h: 0.01
  ny: 19
  nx: 84
  nvar_in: 2
  nvar_out: 1

3.2 Data Reading

The dataset used in this case is stored in the .npz file, and the following code is used to read it.

15
16
17
18
19
20
21
def train(cfg: DictConfig):
    data = np.load(cfg.data_dir)
    coords = data["coords"]
    jinvs = data["jinvs"]
    dxdxis = data["dxdxis"]
    dydxis = data["dydxis"]
    dxdetas = data["dxdetas"]

3.3 Output Transformation Function Construction

This article is a forced boundary constraint. During training, the corresponding output transformation function is used to calculate the differential of the output result of the model.

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
)
sup_constraint = {sup_constraint_res.name: sup_constraint_res}

def _transform_out(
    _input: Dict[str, paddle.Tensor],
    _output: Dict[str, paddle.Tensor],
    pad_singleside: int = cfg.MODEL.pad_singleside,
):
    """Calculation residual.

    Args:
        _input (Dict[str, paddle.Tensor]): The input of the model.
        _output (Dict[str, paddle.Tensor]): The output of the model.
        pad_singleside (int, optional): Pad size. Defaults to cfg.MODEL.pad_singleside.
    """
    output_v = _output["output_v"]
    jinv = _input["jinvs"]
    dxdxi = _input["dxdxis"]
    dydxi = _input["dydxis"]
    dxdeta = _input["dxdetas"]
    dydeta = _input["dydetas"]
    output_v[:, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
    output_v[:, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
    output_v[:, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
    output_v[:, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
    output_v[:, 0, 0, 0] = 0.5 * (output_v[:, 0, 0, 1] + output_v[:, 0, 1, 0])
    output_v[:, 0, 0, -1] = 0.5 * (output_v[:, 0, 0, -2] + output_v[:, 0, 1, -1])
    dvdx = utils.dfdx(output_v, dydeta, dydxi, jinv)
    d2vdx2 = utils.dfdx(dvdx, dydeta, dydxi, jinv)
    dvdy = utils.dfdy(output_v, dxdxi, dxdeta, jinv)

3.4 Constraint Construction

Construct corresponding constraint conditions. Since the boundary constraint is a forced constraint, the constraint conditions are mainly residual constraints.

28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
iters_per_epoch = coords.shape[0]
sup_constraint_res = ppsci.constraint.SupervisedConstraint(
    {
        "dataset": {
            "name": "NamedArrayDataset",
            "input": {
                "coords": coords,
                "jinvs": jinvs,
                "dxdxis": dxdxis,
                "dydxis": dydxis,
                "dxdetas": dxdetas,
                "dydetas": dydetas,
            },
        },
        "batch_size": cfg.TRAIN.batch_size,
        "iters_per_epoch": iters_per_epoch,
        "num_workers": 0,
    },
    ppsci.loss.FunctionalLoss(
        lambda out, label, weight: {"residual": out["residual"]}
    ),

3.5 Optimizer Construction

Consistent with the description in the paper, we use a constant learning rate of 0.001 to construct the Adam optimizer.

25

3.6 Model Training

After completing the above settings, you only need to pass the instantiated objects to ppsci.solver.Solver.

82
83
84
85
86
87
88
89
    return {"residual": paddle.mean(continuity**2)}

model.register_output_transform(_transform_out)
solver = ppsci.solver.Solver(
    model,
    sup_constraint,
    cfg.output_dir,
    optimizer,

Finally start training:

90
epochs=cfg.epochs,

3.7 Model Evaluation

After the model training is completed, the evaluate() function can be used to evaluate and visualize the trained model.

 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    solver.plot_loss_history()


def evaluate(cfg: DictConfig):
    data = np.load(cfg.data_dir)
    coords = data["coords"]

    ofv_sb = paddle.to_tensor(data["OFV_sb"])

    ## create model
    pad_singleside = cfg.MODEL.pad_singleside
    model = ppsci.arch.USCNN(**cfg.MODEL)
    solver = ppsci.solver.Solver(
        model,
        pretrained_model_path=cfg.EVAL.pretrained_model_path,  ### the path of the model
    )
    output_v = solver.predict({"coords": paddle.to_tensor(coords)})
    output_v = output_v["output_v"]

    output_v[0, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
    output_v[0, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
    output_v[0, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
    output_v[0, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
    output_v[0, 0, 0, 0] = 0.5 * (output_v[0, 0, 0, 1] + output_v[0, 0, 1, 0])
    output_v[0, 0, 0, -1] = 0.5 * (output_v[0, 0, 0, -2] + output_v[0, 0, 1, -1])

    ev = paddle.sqrt(
        paddle.mean((ofv_sb - output_v[0, 0]) ** 2) / paddle.mean(ofv_sb**2)
    ).item()
    logger.info(f"ev: {ev}")

    output_v = output_v.numpy()
    ofv_sb = ofv_sb.numpy()
    fig = plt.figure()
    ax = plt.subplot(1, 2, 1)
    utils.visualize(
        ax,
        coords[0, 0, 1:-1, 1:-1],
        coords[0, 1, 1:-1, 1:-1],
        output_v[0, 0, 1:-1, 1:-1],
        "horizontal",
        [0, 1],
    )
    utils.set_axis_label(ax, "p")
    ax.set_title("CNN " + r"$T$")
    ax.set_aspect("equal")
    ax = plt.subplot(1, 2, 2)
    utils.visualize(
        ax,
        coords[0, 0, 1:-1, 1:-1],
        coords[0, 1, 1:-1, 1:-1],
        ofv_sb[1:-1, 1:-1],
        "horizontal",
        [0, 1],
    )
    utils.set_axis_label(ax, "p")
    ax.set_aspect("equal")
    ax.set_title("FV " + r"$T$")

4. Complete Code

heat_equation.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
import os.path as osp
from typing import Dict

import hydra
import matplotlib.pyplot as plt
import numpy as np
import paddle
import utils
from omegaconf import DictConfig

import ppsci
from ppsci.utils import logger


def train(cfg: DictConfig):
    data = np.load(cfg.data_dir)
    coords = data["coords"]
    jinvs = data["jinvs"]
    dxdxis = data["dxdxis"]
    dydxis = data["dydxis"]
    dxdetas = data["dxdetas"]
    dydetas = data["dydetas"]

    model = ppsci.arch.USCNN(**cfg.MODEL)

    optimizer = ppsci.optimizer.Adam(cfg.TRAIN.learning_rate)(model)

    iters_per_epoch = coords.shape[0]
    sup_constraint_res = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": {
                    "coords": coords,
                    "jinvs": jinvs,
                    "dxdxis": dxdxis,
                    "dydxis": dydxis,
                    "dxdetas": dxdetas,
                    "dydetas": dydetas,
                },
            },
            "batch_size": cfg.TRAIN.batch_size,
            "iters_per_epoch": iters_per_epoch,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(
            lambda out, label, weight: {"residual": out["residual"]}
        ),
        name="residual",
    )
    sup_constraint = {sup_constraint_res.name: sup_constraint_res}

    def _transform_out(
        _input: Dict[str, paddle.Tensor],
        _output: Dict[str, paddle.Tensor],
        pad_singleside: int = cfg.MODEL.pad_singleside,
    ):
        """Calculation residual.

        Args:
            _input (Dict[str, paddle.Tensor]): The input of the model.
            _output (Dict[str, paddle.Tensor]): The output of the model.
            pad_singleside (int, optional): Pad size. Defaults to cfg.MODEL.pad_singleside.
        """
        output_v = _output["output_v"]
        jinv = _input["jinvs"]
        dxdxi = _input["dxdxis"]
        dydxi = _input["dydxis"]
        dxdeta = _input["dxdetas"]
        dydeta = _input["dydetas"]
        output_v[:, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
        output_v[:, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
        output_v[:, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
        output_v[:, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
        output_v[:, 0, 0, 0] = 0.5 * (output_v[:, 0, 0, 1] + output_v[:, 0, 1, 0])
        output_v[:, 0, 0, -1] = 0.5 * (output_v[:, 0, 0, -2] + output_v[:, 0, 1, -1])
        dvdx = utils.dfdx(output_v, dydeta, dydxi, jinv)
        d2vdx2 = utils.dfdx(dvdx, dydeta, dydxi, jinv)
        dvdy = utils.dfdy(output_v, dxdxi, dxdeta, jinv)
        d2vdy2 = utils.dfdy(dvdy, dxdxi, dxdeta, jinv)
        continuity = d2vdy2 + d2vdx2
        return {"residual": paddle.mean(continuity**2)}

    model.register_output_transform(_transform_out)
    solver = ppsci.solver.Solver(
        model,
        sup_constraint,
        cfg.output_dir,
        optimizer,
        epochs=cfg.epochs,
        iters_per_epoch=iters_per_epoch,
    )
    solver.train()
    solver.plot_loss_history()


def evaluate(cfg: DictConfig):
    data = np.load(cfg.data_dir)
    coords = data["coords"]

    ofv_sb = paddle.to_tensor(data["OFV_sb"])

    ## create model
    pad_singleside = cfg.MODEL.pad_singleside
    model = ppsci.arch.USCNN(**cfg.MODEL)
    solver = ppsci.solver.Solver(
        model,
        pretrained_model_path=cfg.EVAL.pretrained_model_path,  ### the path of the model
    )
    output_v = solver.predict({"coords": paddle.to_tensor(coords)})
    output_v = output_v["output_v"]

    output_v[0, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
    output_v[0, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
    output_v[0, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
    output_v[0, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
    output_v[0, 0, 0, 0] = 0.5 * (output_v[0, 0, 0, 1] + output_v[0, 0, 1, 0])
    output_v[0, 0, 0, -1] = 0.5 * (output_v[0, 0, 0, -2] + output_v[0, 0, 1, -1])

    ev = paddle.sqrt(
        paddle.mean((ofv_sb - output_v[0, 0]) ** 2) / paddle.mean(ofv_sb**2)
    ).item()
    logger.info(f"ev: {ev}")

    output_v = output_v.numpy()
    ofv_sb = ofv_sb.numpy()
    fig = plt.figure()
    ax = plt.subplot(1, 2, 1)
    utils.visualize(
        ax,
        coords[0, 0, 1:-1, 1:-1],
        coords[0, 1, 1:-1, 1:-1],
        output_v[0, 0, 1:-1, 1:-1],
        "horizontal",
        [0, 1],
    )
    utils.set_axis_label(ax, "p")
    ax.set_title("CNN " + r"$T$")
    ax.set_aspect("equal")
    ax = plt.subplot(1, 2, 2)
    utils.visualize(
        ax,
        coords[0, 0, 1:-1, 1:-1],
        coords[0, 1, 1:-1, 1:-1],
        ofv_sb[1:-1, 1:-1],
        "horizontal",
        [0, 1],
    )
    utils.set_axis_label(ax, "p")
    ax.set_aspect("equal")
    ax.set_title("FV " + r"$T$")
    fig.tight_layout(pad=1)
    fig.savefig(f"{cfg.output_dir}/result.png", bbox_inches="tight")
    plt.close(fig)


def export(cfg: DictConfig):
    model = ppsci.arch.USCNN(**cfg.MODEL)
    # initialize solver
    solver = ppsci.solver.Solver(
        model,
        pretrained_model_path=cfg.INFER.pretrained_model_path,
    )
    # export model
    from paddle.static import InputSpec

    input_spec = [
        {
            key: InputSpec([None, 2, 19, 84], "float32", name=key)
            for key in model.input_keys
        },
    ]
    solver.export(input_spec, cfg.INFER.export_path)


def inference(cfg: DictConfig):
    from deploy.python_infer import pinn_predictor

    predictor = pinn_predictor.PINNPredictor(cfg)
    data = np.load(cfg.data_dir)
    coords = data["coords"]
    ofv_sb = data["OFV_sb"]

    ## create model
    pad_singleside = cfg.MODEL.pad_singleside
    input_spec = {"coords": coords}

    output_v = predictor.predict(input_spec, cfg.INFER.batch_size)
    # mapping data to cfg.INFER.output_keys
    output_v = {
        store_key: output_v[infer_key]
        for store_key, infer_key in zip(cfg.MODEL.output_keys, output_v.keys())
    }

    output_v = output_v["output_v"]

    output_v[0, 0, -pad_singleside:, pad_singleside:-pad_singleside] = 0
    output_v[0, 0, :pad_singleside, pad_singleside:-pad_singleside] = 1
    output_v[0, 0, pad_singleside:-pad_singleside, -pad_singleside:] = 1
    output_v[0, 0, pad_singleside:-pad_singleside, 0:pad_singleside] = 1
    output_v[0, 0, 0, 0] = 0.5 * (output_v[0, 0, 0, 1] + output_v[0, 0, 1, 0])
    output_v[0, 0, 0, -1] = 0.5 * (output_v[0, 0, 0, -2] + output_v[0, 0, 1, -1])

    ev = np.sqrt(np.mean((ofv_sb - output_v[0, 0]) ** 2) / np.mean(ofv_sb**2))
    logger.info(f"ev: {ev}")

    fig = plt.figure()
    ax = plt.subplot(1, 2, 1)
    utils.visualize(
        ax,
        coords[0, 0, 1:-1, 1:-1],
        coords[0, 1, 1:-1, 1:-1],
        output_v[0, 0, 1:-1, 1:-1],
        "horizontal",
        [0, 1],
    )
    utils.set_axis_label(ax, "p")
    ax.set_title("CNN " + r"$T$")
    ax.set_aspect("equal")
    ax = plt.subplot(1, 2, 2)
    utils.visualize(
        ax,
        coords[0, 0, 1:-1, 1:-1],
        coords[0, 1, 1:-1, 1:-1],
        ofv_sb[1:-1, 1:-1],
        "horizontal",
        [0, 1],
    )
    utils.set_axis_label(ax, "p")
    ax.set_aspect("equal")
    ax.set_title("FV " + r"$T$")
    fig.tight_layout(pad=1)
    fig.savefig(osp.join(cfg.output_dir, "result.png"), bbox_inches="tight")
    plt.close(fig)


@hydra.main(version_base=None, config_path="./conf", config_name="heat_equation.yaml")
def main(cfg: DictConfig):
    if cfg.mode == "train":
        train(cfg)
    elif cfg.mode == "eval":
        evaluate(cfg)
    elif cfg.mode == "export":
        export(cfg)
    elif cfg.mode == "infer":
        inference(cfg)
    else:
        raise ValueError(
            f"cfg.mode should in ['train', 'eval', 'export', 'infer'], but got '{cfg.mode}'"
        )


if __name__ == "__main__":
    main()

5. Result Display

Heat equation result display: image

Heat equation with boundary result display:

T=0 image

T=3 image

T=6 image

6. Summary

This paper constructs a coordinate transformation function using harmonic mapping, so that the physical information network can be trained on irregular non-uniform grids. At the same time, because the mapping is performed using traditional methods, it can be embedded before and after the network without training. Through a large number of experiments, it is shown that the network can perform better than SOAT networks on various irregular grid problems.

7. References

PhyGeoNet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state PDEs on irregular domain

Github PhyGeoNet