跳转至

Euler Beam

python euler_beam.py
python euler_beam.py mode=eval EVAL.pretrained_model_path=https://paddle-org.bj.bcebos.com/paddlescience/models/euler_beam/euler_beam_pretrained.pdparams
python euler_beam.py mode=export
python euler_beam.py mode=infer
预训练模型 指标
euler_beam_pretrained.pdparams loss(L2Rel_Metric): 0.00000
L2Rel.u(L2Rel_Metric): 0.00058

1. 问题定义

Euler Beam 公式:

\[ \dfrac{\partial^{4} u}{\partial x^{4}} + 1 = 0, x \in [0, 1] \]

边界条件:

\[ u''(1)=0, u'''(1)=0 \]

狄利克雷条件:

\[ u(0)=0 \]

诺依曼边界条件:

\[ u'(0)=0 \]

2. 问题求解

接下来开始讲解如何将问题一步一步地转化为 PaddleScience 代码,用深度学习的方法求解该问题。 为了快速理解 PaddleScience,接下来仅对模型构建、方程构建、计算域构建等关键步骤进行阐述,而其余细节请参考 API文档

2.1 模型构建

在 Euler Beam 问题中,每一个已知的坐标点 \(x\) 都有对应的待求解的未知量 \(u\) ,我们在这里使用比较简单的 MLP(Multilayer Perceptron, 多层感知机) 来表示 \(x\)\(u\) 的映射函数 \(f: \mathbb{R}^1 \to \mathbb{R}^1\) ,即:

\[ u = f(x) \]

上式中 \(f\) 即为 MLP 模型本身,用 PaddleScience 代码表示如下

# set model
model = ppsci.arch.MLP(**cfg.MODEL)

其中,用于初始化模型的参数通过配置文件进行配置:

# model settings
MODEL:
  input_keys: ["x"]
  output_keys: ["u"]
  num_layers: 3
  hidden_size: 20

接着通过指定 MLP 的层数、神经元个数,我们就实例化出了一个拥有 3 层隐藏神经元,每层神经元数为 20 的神经网络模型 model

2.2 方程构建

Euler Beam 的方程构建可以直接使用 PaddleScience 内置的 Biharmonic,指定该类的参数 dim 为 1,q 为 -1,D 为1。

# set equation(s)
equation = {"biharmonic": ppsci.equation.Biharmonic(dim=1, q=cfg.q, D=cfg.D)}

2.3 计算域构建

本文中 Euler Beam 问题作用在以 (0.0, 1.0) 的一维区域上, 因此可以直接使用 PaddleScience 内置的空间几何 Interval 作为计算域。

# set geometry
geom = {"interval": ppsci.geometry.Interval(0, 1)}

2.4 约束构建

在本案例中,我们使用了两个约束条件在计算域中指导模型的训练分别是作用于采样点上的方程约束和作用于边界点上的约束。

在定义约束之前,需要给每一种约束指定采样点个数,表示每一种约束在其对应计算域内采样数据的数量,以及通用的采样配置。

# training settings
TRAIN:
  epochs: 10000
  iters_per_epoch: 1
  save_freq: 1000
  eval_during_train: true
  eval_freq: 1000
  learning_rate: 1.0e-3
  batch_size:
    pde: 100
    bc: 4

2.4.1 内部点约束

以作用在内部点上的 InteriorConstraint 为例,代码如下:

# set dataloader config
dataloader_cfg = {
    "dataset": "IterableNamedArrayDataset",
    "iters_per_epoch": cfg.TRAIN.iters_per_epoch,
}
# set constraint
pde_constraint = ppsci.constraint.InteriorConstraint(
    equation["biharmonic"].equations,
    {"biharmonic": 0},
    geom["interval"],
    {**dataloader_cfg, "batch_size": cfg.TRAIN.batch_size.pde},
    ppsci.loss.MSELoss(),
    random="Hammersley",
    name="EQ",
)

2.4.2 边界约束

同理,我们还需要构建边界的约束。但与构建 InteriorConstraint 约束不同的是,由于作用区域是边界,因此我们使用 BoundaryConstraint 类,代码如下:

bc = ppsci.constraint.BoundaryConstraint(
    {
        "u0": lambda d: d["u"][0:1],
        "u__x": lambda d: jacobian(d["u"], d["x"])[1:2],
        "u__x__x": lambda d: hessian(d["u"], d["x"])[2:3],
        "u__x__x__x": lambda d: jacobian(hessian(d["u"], d["x"]), d["x"])[3:4],
    },
    {"u0": 0, "u__x": 0, "u__x__x": 0, "u__x__x__x": 0},
    geom["interval"],
    {**dataloader_cfg, "batch_size": cfg.TRAIN.batch_size.bc},
    ppsci.loss.MSELoss("sum"),
    evenly=True,
    name="BC",
)

2.5 超参数设定

接下来我们需要在配置文件中指定训练轮数,此处我们按实验经验,使用一万轮训练轮数,评估间隔为一千轮。

# training settings
TRAIN:
  epochs: 10000
  iters_per_epoch: 1
  save_freq: 1000
  eval_during_train: true
  eval_freq: 1000

2.6 优化器构建

训练过程会调用优化器来更新模型参数,此处选择较为常用的 Adam 优化器。

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

2.7 评估器构建

在训练过程中通常会按一定轮数间隔,用验证集(测试集)评估当前模型的训练情况,因此使用 ppsci.validate.GeometryValidator 构建评估器。

l2_rel_metric = ppsci.validate.GeometryValidator(
    {"u": lambda out: out["u"]},
    {"u": u_solution_func},
    geom["interval"],
    {
        "dataset": "IterableNamedArrayDataset",
        "total_size": cfg.EVAL.total_size,
    },
    ppsci.loss.MSELoss(),
    evenly=True,
    metric={"L2Rel": ppsci.metric.L2Rel()},
    name="L2Rel_Metric",
)
validator = {l2_rel_metric.name: l2_rel_metric}

2.8 可视化器构建

在模型评估时,如果评估结果是可以可视化的数据,我们可以选择合适的可视化器来对输出结果进行可视化。

本文中的输出数据是一个曲线图,因此我们只需要将评估的输出数据保存成 png 文件即可。代码如下:

# set visualizer(optional)
visu_points = geom["interval"].sample_interior(cfg.EVAL.total_size, evenly=True)
visualizer = {
    "visualize_u": ppsci.visualize.VisualizerScatter1D(
        visu_points,
        ("x",),
        {
            "u_label": lambda d: u_solution_func(d),
            "u_pred": lambda d: d["u"],
        },
        num_timestamps=1,
        prefix="result_u",
    )
}

2.9 模型训练、评估与可视化

完成上述设置之后,只需要将上述实例化的对象按顺序传递给 ppsci.solver.Solver,然后启动训练、评估、可视化。

# initialize solver
solver = ppsci.solver.Solver(
    model,
    constraint,
    cfg.output_dir,
    optimizer,
    epochs=cfg.TRAIN.epochs,
    iters_per_epoch=cfg.TRAIN.iters_per_epoch,
    eval_during_train=cfg.TRAIN.eval_during_train,
    eval_freq=cfg.TRAIN.eval_freq,
    seed=cfg.seed,
    equation=equation,
    geom=geom,
    validator=validator,
    visualizer=visualizer,
    pretrained_model_path=cfg.TRAIN.pretrained_model_path,
    checkpoint_path=cfg.TRAIN.checkpoint_path,
    eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
    to_static=cfg.to_static,
)
# train model
solver.train()
# evaluate after finished training
solver.eval()
# visualize prediction after finished training
solver.visualize()

3. 完整代码

euler_beam.py
# Copyright (c) 2023 PaddlePaddle Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import hydra
from omegaconf import DictConfig

import ppsci
from ppsci.autodiff import hessian
from ppsci.autodiff import jacobian


def train(cfg: DictConfig):
    # set model
    model = ppsci.arch.MLP(**cfg.MODEL)

    # set geometry
    geom = {"interval": ppsci.geometry.Interval(0, 1)}

    # set equation(s)
    equation = {"biharmonic": ppsci.equation.Biharmonic(dim=1, q=cfg.q, D=cfg.D)}

    # set dataloader config
    dataloader_cfg = {
        "dataset": "IterableNamedArrayDataset",
        "iters_per_epoch": cfg.TRAIN.iters_per_epoch,
    }
    # set constraint
    pde_constraint = ppsci.constraint.InteriorConstraint(
        equation["biharmonic"].equations,
        {"biharmonic": 0},
        geom["interval"],
        {**dataloader_cfg, "batch_size": cfg.TRAIN.batch_size.pde},
        ppsci.loss.MSELoss(),
        random="Hammersley",
        name="EQ",
    )
    bc = ppsci.constraint.BoundaryConstraint(
        {
            "u0": lambda d: d["u"][0:1],
            "u__x": lambda d: jacobian(d["u"], d["x"])[1:2],
            "u__x__x": lambda d: hessian(d["u"], d["x"])[2:3],
            "u__x__x__x": lambda d: jacobian(hessian(d["u"], d["x"]), d["x"])[3:4],
        },
        {"u0": 0, "u__x": 0, "u__x__x": 0, "u__x__x__x": 0},
        geom["interval"],
        {**dataloader_cfg, "batch_size": cfg.TRAIN.batch_size.bc},
        ppsci.loss.MSELoss("sum"),
        evenly=True,
        name="BC",
    )
    # wrap constraints together
    constraint = {
        pde_constraint.name: pde_constraint,
        bc.name: bc,
    }

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

    # set validator
    def u_solution_func(out):
        """compute ground truth for u as label data"""
        x = out["x"]
        return -(x**4) / 24 + x**3 / 6 - x**2 / 4

    l2_rel_metric = ppsci.validate.GeometryValidator(
        {"u": lambda out: out["u"]},
        {"u": u_solution_func},
        geom["interval"],
        {
            "dataset": "IterableNamedArrayDataset",
            "total_size": cfg.EVAL.total_size,
        },
        ppsci.loss.MSELoss(),
        evenly=True,
        metric={"L2Rel": ppsci.metric.L2Rel()},
        name="L2Rel_Metric",
    )
    validator = {l2_rel_metric.name: l2_rel_metric}

    # set visualizer(optional)
    visu_points = geom["interval"].sample_interior(cfg.EVAL.total_size, evenly=True)
    visualizer = {
        "visualize_u": ppsci.visualize.VisualizerScatter1D(
            visu_points,
            ("x",),
            {
                "u_label": lambda d: u_solution_func(d),
                "u_pred": lambda d: d["u"],
            },
            num_timestamps=1,
            prefix="result_u",
        )
    }

    # initialize solver
    solver = ppsci.solver.Solver(
        model,
        constraint,
        cfg.output_dir,
        optimizer,
        epochs=cfg.TRAIN.epochs,
        iters_per_epoch=cfg.TRAIN.iters_per_epoch,
        eval_during_train=cfg.TRAIN.eval_during_train,
        eval_freq=cfg.TRAIN.eval_freq,
        seed=cfg.seed,
        equation=equation,
        geom=geom,
        validator=validator,
        visualizer=visualizer,
        pretrained_model_path=cfg.TRAIN.pretrained_model_path,
        checkpoint_path=cfg.TRAIN.checkpoint_path,
        eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
        to_static=cfg.to_static,
    )
    # train model
    solver.train()
    # evaluate after finished training
    solver.eval()
    # visualize prediction after finished training
    solver.visualize()


def evaluate(cfg: DictConfig):
    # set model
    model = ppsci.arch.MLP(**cfg.MODEL)

    # set geometry
    geom = {"interval": ppsci.geometry.Interval(0, 1)}

    # set equation(s)
    equation = {"biharmonic": ppsci.equation.Biharmonic(dim=1, q=cfg.q, D=cfg.D)}

    # set validator
    def u_solution_func(out):
        """compute ground truth for u as label data"""
        x = out["x"]
        return -(x**4) / 24 + x**3 / 6 - x**2 / 4

    l2_rel_metric = ppsci.validate.GeometryValidator(
        {"u": lambda out: out["u"]},
        {"u": u_solution_func},
        geom["interval"],
        {
            "dataset": "IterableNamedArrayDataset",
            "total_size": cfg.EVAL.total_size,
        },
        ppsci.loss.MSELoss(),
        evenly=True,
        metric={"L2Rel": ppsci.metric.L2Rel()},
        name="L2Rel_Metric",
    )
    validator = {l2_rel_metric.name: l2_rel_metric}

    # set visualizer(optional)
    visu_points = geom["interval"].sample_interior(cfg.EVAL.total_size, evenly=True)
    visualizer = {
        "visualize_u": ppsci.visualize.VisualizerScatter1D(
            visu_points,
            ("x",),
            {
                "u_label": lambda d: u_solution_func(d),
                "u_pred": lambda d: d["u"],
            },
            num_timestamps=1,
            prefix="result_u",
        )
    }

    # initialize solver
    solver = ppsci.solver.Solver(
        model,
        None,
        cfg.output_dir,
        None,
        seed=cfg.seed,
        equation=equation,
        geom=geom,
        validator=validator,
        visualizer=visualizer,
        pretrained_model_path=cfg.EVAL.pretrained_model_path,
        eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
        to_static=cfg.to_static,
    )
    # evaluate after finished training
    solver.eval()
    # visualize prediction after finished training
    solver.visualize()


def export(cfg: DictConfig):
    # set model
    model = ppsci.arch.MLP(**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, 1], "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)

    # set geometry
    geom = {"interval": ppsci.geometry.Interval(0, 1)}
    input_dict = geom["interval"].sample_interior(cfg.INFER.total_size, evenly=True)

    output_dict = predictor.predict({"x": input_dict["x"]}, cfg.INFER.batch_size)

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

    def u_solution_func(out):
        """compute ground truth for u as label data"""
        x = out["x"]
        return -(x**4) / 24 + x**3 / 6 - x**2 / 4

    ppsci.visualize.save_plot_from_1d_dict(
        "./euler_beam_pred",
        {**input_dict, **output_dict, "u_label": u_solution_func(input_dict)},
        ("x",),
        ("u", "u_label"),
    )


@hydra.main(version_base=None, config_path="./conf", config_name="euler_beam.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()

4. 结果展示

使用训练得到的模型对上述计算域中均匀取的共 NPOINT_TOTAL 个点 \(x_i\) 进行预测,预测结果如下所示。图像中横坐标为 \(x\),纵坐标为对应的预测结果 \(u\)

euler_beam

模型预测结果