Euler Beam

Model Training CommandModel Evaluation CommandModel Export CommandModel Inference Command

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

Pretrained Model	Metrics
euler_beam_pretrained.pdparams	loss(L2Rel_Metric): 0.00000 L2Rel.u(L2Rel_Metric): 0.00058

1. Problem Definition

Euler Beam Formula:

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

Boundary Conditions:

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

Dirichlet Condition:

\[ u(0)=0 \]

Neumann Boundary Condition:

\[ u'(0)=0 \]

2. Problem Solving

Next, we will explain how to convert the problem into PaddleScience code step by step and solve the problem using deep learning methods. In order to quickly understand PaddleScience, only key steps such as model construction, equation construction, and computational domain construction are described below, while other details please refer to API Documentation.

2.1 Model Construction

In the Euler Beam problem, each known coordinate point \(x\) has a corresponding unknown quantity \(u\) to be solved. We use a relatively simple MLP (Multilayer Perceptron) here to represent the mapping function \(f: \mathbb{R}^1 \to \mathbb{R}^1\) from \(x\) to \(u\), i.e.:

\[ u = f(x) \]

In the above formula, \(f\) is the MLP model itself, expressed in PaddleScience code as follows:

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

The parameters used to initialize the model are configured through the configuration file:

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

Then by specifying the number of layers and neurons of MLP, we instantiated a neural network model model with 3 hidden layers and 20 neurons per layer.

2.2 Equation Construction

The equation construction of Euler Beam can directly use the Biharmonic built in PaddleScience, specifying the parameter dim of this class as 1, q as -1, and D as 1.

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

2.3 Computational Domain Construction

In this article, the Euler Beam problem acts on a one-dimensional area of (0.0, 1.0), so the spatial geometry Interval built in PaddleScience can be used directly as the computational domain.

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

2.4 Constraint Construction

In this case, we used two constraints to guide the training of the model in the computational domain, namely the equation constraint acting on the sampling point and the constraint acting on the boundary point.

Before defining constraints, you need to specify the number of sampling points for each constraint, indicating the number of sampled data for each constraint in its corresponding computational domain, as well as general sampling configuration.

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
  pretrained_model_path: null

2.4.1 Interior Point Constraint

Taking InteriorConstraint acting on internal points as an example, the code is as follows:

# 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 Boundary Constraint

Similarly, we also need to construct boundary constraints. However, unlike constructing InteriorConstraint, since the action area is the boundary, we use the BoundaryConstraint class, code as follows:

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 Hyperparameter Setting

Next, we need to specify the number of training epochs in the configuration file. Here, based on experimental experience, we use 10,000 training epochs, with an evaluation interval of 1,000 epochs.

TRAIN:
  epochs: 10000
  iters_per_epoch: 1
  save_freq: 1000
  eval_during_train: true
  eval_freq: 1000
  learning_rate: 1.0e-3

2.6 Optimizer Construction

The training process will call the optimizer to update model parameters. Here, the more commonly used Adam optimizer is selected.

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

2.7 Validator Construction

Usually during the training process, the training status of the current model is evaluated using the validation set (test set) at a certain epoch interval, so ppsci.validate.GeometryValidator is used to construct the validator.

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 Visualizer Construction

During model evaluation, if the evaluation result is data that can be visualized, we can choose a suitable visualizer to visualize the output result.

The output data in this article is a graph, so we only need to save the evaluated output data as a png file. The code is as follows:

# 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 Model Training, Evaluation and Visualization

After completing the above settings, you only need to pass the instantiated objects to ppsci.solver.Solver in order, and then start training, evaluation, and visualization.

# 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. Complete Code

# 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. Result Display

The trained model is used to predict a total of NPOINT_TOTAL points \(x_i\) uniformly taken from the above computational domain. The prediction results are shown below. In the image, the abscissa is \(x\), and the ordinate is the corresponding predicted result \(u\).

Model prediction result