Skip to content

2D-Cylinder (2D Flow Around a Cylinder)

AI Studio Quick Experience

# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar -o cylinder2d_unsteady_Re100_dataset.tar
# unzip it
tar -xvf cylinder2d_unsteady_Re100_dataset.tar
python cylinder2d_unsteady_Re100.py

# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar -o cylinder2d_unsteady_Re100_dataset.tar
# unzip it
tar -xvf cylinder2d_unsteady_Re100_dataset.tar
python cylinder2d_unsteady_Re100.py mode=eval EVAL.pretrained_model_path=https://paddle-org.bj.bcebos.com/paddlescience/models/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_pretrained.pdparams

python cylinder2d_unsteady_Re100.py mode=export

# linux
wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar -o cylinder2d_unsteady_Re100_dataset.tar
# unzip it
tar -xvf cylinder2d_unsteady_Re100_dataset.tar
python cylinder2d_unsteady_Re100.py mode=infer
Pretrained Model Metrics
cylinder2d_unsteady_Re100_pretrained.pdparams loss(Residual): 0.00398
MSE.continuity(Residual): 0.00126
MSE.momentum_x(Residual): 0.00151
MSE.momentum_y(Residual): 0.00120

1. Background Introduction

The flow around a cylinder is a fundamental problem with applications in diverse fields. In industrial design, it is used to simulate and optimize fluid dynamics in equipment such as wind turbines, as well as the aerodynamic performance of automobiles and aircraft. In environmental engineering, it aids in predicting river flooding and modeling pollutant dispersion. Furthermore, it holds significant practical value in broader engineering contexts, including fluid dynamics, heat transfer, and aerodynamics.

2D flow around a cylinder describes the flow pattern of low-speed steady flow past a two-dimensional cylinder, governed primarily by the Reynolds number (\(Re\)). When \(Re \le 1\) (the Stokes region), inertial forces are negligible compared to viscous forces; streamlines are symmetrical upstream and downstream, and the drag coefficient is approximately inversely proportional to \(Re\) (ranging from 10 to 60). As \(Re\) increases, this symmetry is gradually lost.

2. Problem Definition

Mass conservation:

\[ \dfrac{\partial u}{\partial x} + \dfrac{\partial v}{\partial y} = 0 \]

\(x\) momentum conservation:

\[ \dfrac{\partial u}{\partial t} + u\dfrac{\partial u}{\partial x} + v\dfrac{\partial u}{\partial y} = -\dfrac{1}{\rho}\dfrac{\partial p}{\partial x} + \nu(\dfrac{\partial ^2 u}{\partial x ^2} + \dfrac{\partial ^2 u}{\partial y ^2}) \]

\(y\) momentum conservation:

\[ \dfrac{\partial v}{\partial t} + u\dfrac{\partial v}{\partial x} + v\dfrac{\partial v}{\partial y} = -\dfrac{1}{\rho}\dfrac{\partial p}{\partial y} + \nu(\dfrac{\partial ^2 v}{\partial x ^2} + \dfrac{\partial ^2 v}{\partial y ^2}) \]

Let:

\(t^* = \dfrac{L}{U_0}\)

\(x^*=y^* = L\)

\(u^*=v^* = U_0\)

\(p^* = \rho {U_0}^2\)

Define:

Dimensionless time \(\tau = \dfrac{t}{t^*}\)

Dimensionless coordinate \(x: X = \dfrac{x}{x^*}\); Dimensionless coordinate \(y: Y = \dfrac{y}{y^*}\)

Dimensionless velocity \(x: U = \dfrac{u}{u^*}\); Dimensionless velocity \(y: V = \dfrac{v}{u^*}\)

Dimensionless pressure \(P = \dfrac{p}{p^*}\)

Reynolds number \(Re = \dfrac{L U_0}{\nu}\)

The following dimensionless Navier-Stokes equations apply to the interior of the fluid domain:

Mass conservation:

\[ \dfrac{\partial U}{\partial X} + \dfrac{\partial U}{\partial Y} = 0 \]

\(x\) momentum conservation:

\[ \dfrac{\partial U}{\partial \tau} + U\dfrac{\partial U}{\partial X} + V\dfrac{\partial U}{\partial Y} = -\dfrac{\partial P}{\partial X} + \dfrac{1}{Re}(\dfrac{\partial ^2 U}{\partial X^2} + \dfrac{\partial ^2 U}{\partial Y^2}) \]

\(y\) momentum conservation:

\[ \dfrac{\partial V}{\partial \tau} + U\dfrac{\partial V}{\partial X} + V\dfrac{\partial V}{\partial Y} = -\dfrac{\partial P}{\partial Y} + \dfrac{1}{Re}(\dfrac{\partial ^2 V}{\partial X^2} + \dfrac{\partial ^2 V}{\partial Y^2}) \]

For the fluid domain boundary and the inner circumference boundary of the fluid domain, Dirichlet boundary conditions need to be applied:

Fluid domain inlet boundary:

\[ u=1, v=0 \]

Circumference boundary:

\[ u=0, v=0 \]

Fluid domain outlet boundary:

\[ p=0 \]

3. 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.

Before starting to build the code, please download the dataset required for training and evaluation according to the following command

wget -c https://paddle-org.bj.bcebos.com/paddlescience/datasets/cylinder2d_unsteady_Re100/cylinder2d_unsteady_Re100_dataset.tar
tar -xf cylinder2d_unsteady_Re100_dataset.tar

3.1 Model Construction

In the 2D-Cylinder problem, for each spatiotemporal coordinate \((t, x, y)\), there are three unknown quantities to solve: lateral velocity \(u\), longitudinal velocity \(v\), and pressure \(p\). We employ a Multilayer Perceptron (MLP) to approximate the mapping function \(f: \mathbb{R}^3 \to \mathbb{R}^3\) from \((t, x, y)\) to \((u, v, p)\), such that:

\[ u, v, p = f(t, x, y) \]

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

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

To ensure accurate and efficient variable access during computation, we define the model's input keys as ["t", "x", "y"] and output keys as ["u", "v", "p"], maintaining consistency with the code.

We instantiate the MLP model with 9 hidden layers, 50 neurons per layer, and the tanh activation function.

3.2 Equation Construction

Since this problem involves the 2D transient Navier-Stokes equations, we can directly use the built-in NavierStokes class in PaddleScience.

36
37
38
39
# set equation
equation = {
    "NavierStokes": ppsci.equation.NavierStokes(cfg.VISCOSITY, cfg.DENSITY, 2, True)
}

When instantiating the NavierStokes class, necessary parameters need to be specified: dynamic viscosity \(\nu=0.02\), fluid density \(\rho=1.0\).

3.3 Computational Domain Construction

The computational domain is defined by point clouds stored in CSV files. We combine the built-in PointCloud geometry and TimeDomain to construct a spatiotemporal TimeXGeometry domain.

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
# set timestamps
train_timestamps = np.linspace(
    cfg.TIME_START, cfg.TIME_END, cfg.NUM_TIMESTAMPS, endpoint=True
).astype("float32")
train_timestamps = np.random.choice(train_timestamps, cfg.TRAIN_NUM_TIMESTAMPS)
train_timestamps.sort()
t0 = np.array([cfg.TIME_START], dtype="float32")

val_timestamps = np.linspace(
    cfg.TIME_START, cfg.TIME_END, cfg.NUM_TIMESTAMPS, endpoint=True
).astype("float32")

logger.message(f"train_timestamps: {train_timestamps.tolist()}")
logger.message(f"val_timestamps: {val_timestamps.tolist()}")

# set time-geometry
geom = {
    "time_rect": ppsci.geometry.TimeXGeometry(
        ppsci.geometry.TimeDomain(
            cfg.TIME_START,
            cfg.TIME_END,
            timestamps=np.concatenate((t0, train_timestamps), axis=0),
        ),
        ppsci.geometry.PointCloud(
            reader.load_csv_file(
                "./datasets/domain_train.csv",
                ("x", "y"),
                alias_dict={"x": "Points:0", "y": "Points:1"},
            ),
            ("x", "y"),
        ),
    ),
    "time_rect_eval": ppsci.geometry.PointCloud(
        reader.load_csv_file(
            "./datasets/domain_eval.csv",
            ("t", "x", "y"),
        ),
        ("t", "x", "y"),
    ),
}
  1. The evaluation data points already contain timestamp information, so there is no need to combine with TimeDomain into TimeXGeometry, just use PointCloud to read in the data.
Tip

PointCloud and TimeDomain are two Geometry derived classes that can be used independently.

If the input data only comes from point cloud geometry, you can directly use ppsci.geometry.PointCloud(...) to create a spatial geometric domain object;

If the input data only comes from a one-dimensional time domain, you can directly use ppsci.geometry.TimeDomain(...) to construct a time domain object.

3.4 Constraint Construction

According to the dimensionless formulas and boundary conditions obtained in 2. Problem Definition, corresponding to the three constraint conditions guiding model training in the computational domain, namely:

  1. Dimensionless Navier-Stokes equation constraint applied to internal points of the fluid domain (after simple term shifting)

    \[ \dfrac{\partial U}{\partial X} + \dfrac{\partial U}{\partial Y} = 0 \]
    \[ \dfrac{\partial U}{\partial \tau} + U\dfrac{\partial U}{\partial X} + V\dfrac{\partial U}{\partial Y} + \dfrac{\partial P}{\partial X} - \dfrac{1}{Re}(\dfrac{\partial ^2 U}{\partial X^2} + \dfrac{\partial ^2 U}{\partial Y^2}) = 0 \]
    \[ \dfrac{\partial V}{\partial \tau} + U\dfrac{\partial V}{\partial X} + V\dfrac{\partial V}{\partial Y} + \dfrac{\partial P}{\partial Y} - \dfrac{1}{Re}(\dfrac{\partial ^2 V}{\partial X^2} + \dfrac{\partial ^2 V}{\partial Y^2}) = 0 \]

    In order to facilitate obtaining intermediate variables, the NavierStokes class internally names the results on the left side of the above formula as continuity, momentum_x, momentum_y respectively.

  2. Dirichlet boundary condition constraints applied to the fluid domain inlet, internal circumference, and fluid domain outlet

    Fluid domain inlet boundary:

    \[ u=1, v=0 \]

    Fluid domain outlet boundary:

    \[ p=0 \]

    Circumference boundary:

    \[ u=0, v=0 \]

  3. Initial value condition constraint applied to internal points of the fluid domain at the initial moment:

    \[ u=u_{t0}, v=v_{t0}, p=p_{t0} \]

Next, use InteriorConstraint and SupervisedConstraint built in PaddleScience to construct the above two constraints.

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.

82
83
84
# pde/bc/sup constraint use t1~tn, initial constraint use t0
NTIME_PDE = len(train_timestamps)
ALIAS_DICT = {"x": "Points:0", "y": "Points:1", "u": "U:0", "v": "U:1"}

3.4.1 Interior Point Constraint

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

86
87
88
89
90
91
92
93
94
95
96
97
98
# set constraint
pde_constraint = ppsci.constraint.InteriorConstraint(
    equation["NavierStokes"].equations,
    {"continuity": 0, "momentum_x": 0, "momentum_y": 0},
    geom["time_rect"],
    {
        "dataset": "IterableNamedArrayDataset",
        "batch_size": cfg.NPOINT_PDE * NTIME_PDE,
        "iters_per_epoch": cfg.TRAIN.iters_per_epoch,
    },
    ppsci.loss.MSELoss("mean"),
    name="EQ",
)
  • Equation Expression: Specifies how to calculate the constraint target. We use equation["NavierStokes"].equations from Section 3.2.
  • Target Values: The target values for the constraint variables. We aim to minimize the residuals of continuity, momentum_x, and momentum_y to 0.
  • Computational Domain: The domain where the constraint applies. We use geom["time_rect"] from Section 3.3.
  • Sampling Configuration: Specifies how data is sampled. We use full data points (IterableNamedArrayDataset, iters_per_epoch=1) with a batch_size of 9420 * 30 (9420 spatial points × 30 timestamps).
  • Loss Function: We use the Mean Squared Error (MSE) with reduction="mean" to average the loss across all points.
  • Name: A unique name for the constraint, set here as "EQ".

3.4.2 Boundary Constraint

We also construct Dirichlet boundary constraints for the inflow, outflow, and cylinder circumference boundaries. Taking bc_inlet_cylinder as an example, we use SupervisedConstraint since the boundary data is provided in a CSV file. The dataloader_cfg is configured as follows:

  • Dataset Class: IterableCSVDataset for loading full data.
  • File Path: ./datasets/domain_inlet_cylinder.csv.
  • Input Keys: ("x", "y") to be read from the file.
  • Label Keys: ("u", "v") to be read from the file.
  • Alias Dict: Maps file column names to standard keys, e.g., {"x": "Points:0", "y": "Points:1", "u": "U:0", "v": "U:1"}.
  • Weight Dict: Assigns weights to labels. We amplify "u" and "v" weights to 10: {"u": 10, "v": 10}.
  • Timestamps: Specifies the time information, set here to train_timestamps.

We use the MSE loss function with reduction="mean" and name the constraint "BC_inlet_cylinder".

The remaining bc_outlet is constructed according to the same principle, the code is as follows:

 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
bc_inlet_cylinder = ppsci.constraint.SupervisedConstraint(
    {
        "dataset": {
            "name": "IterableCSVDataset",
            "file_path": cfg.DOMAIN_INLET_CYLINDER_PATH,
            "input_keys": ("x", "y"),
            "label_keys": ("u", "v"),
            "alias_dict": ALIAS_DICT,
            "weight_dict": {"u": 10, "v": 10},
            "timestamps": train_timestamps,
        },
    },
    ppsci.loss.MSELoss("mean"),
    name="BC_inlet_cylinder",
)
bc_outlet = ppsci.constraint.SupervisedConstraint(
    {
        "dataset": {
            "name": "IterableCSVDataset",
            "file_path": cfg.DOMAIN_OUTLET_PATH,
            "input_keys": ("x", "y"),
            "label_keys": ("p",),
            "alias_dict": ALIAS_DICT,
            "timestamps": train_timestamps,
        },
    },
    ppsci.loss.MSELoss("mean"),
    name="BC_outlet",
)

3.4.3 Initial Value Constraint

For points in the fluid domain at time \(t=t_0\), we also need to apply initial value constraints to \(u\), \(v\), \(p\). The code is as follows:

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
ic = ppsci.constraint.SupervisedConstraint(
    {
        "dataset": {
            "name": "IterableCSVDataset",
            "file_path": cfg.IC0_1_PATH,
            "input_keys": ("x", "y"),
            "label_keys": ("u", "v", "p"),
            "alias_dict": ALIAS_DICT,
            "weight_dict": {"u": 10, "v": 10, "p": 10},
            "timestamps": t0,
        },
    },
    ppsci.loss.MSELoss("mean"),
    name="IC",
)

3.4.4 Supervised Constraint

In this case, a certain number of supervision points are added inside the fluid domain to ensure the final convergence of the model, so a supervised constraint needs to be added finally. The data also comes from CSV files. The code is as follows:

143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
sup_constraint = ppsci.constraint.SupervisedConstraint(
    {
        "dataset": {
            "name": "IterableCSVDataset",
            "file_path": cfg.PROBE1_50_PATH,
            "input_keys": ("t", "x", "y"),
            "label_keys": ("u", "v"),
            "alias_dict": ALIAS_DICT,
            "weight_dict": {"u": 10, "v": 10},
            "timestamps": train_timestamps,
        },
    },
    ppsci.loss.MSELoss("mean"),
    name="Sup",
)

After the differential equation constraint, boundary constraint, initial value constraint, and supervised constraint are constructed, encapsulate them into a dictionary with the names we just named as keys for subsequent access.

159
160
161
162
163
164
165
166
# wrap constraints together
constraint = {
    pde_constraint.name: pde_constraint,
    bc_inlet_cylinder.name: bc_inlet_cylinder,
    bc_outlet.name: bc_outlet,
    ic.name: ic,
    sup_constraint.name: sup_constraint,
}

3.5 Hyperparameter Setting

Next, we need to specify the number of training epochs and learning rate. Here, based on experimental experience, we use 40,000 training epochs, evaluation interval is 400 epochs, and learning rate is set to 0.001.

60

3.6 Optimizer Construction

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

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

3.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.

171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
# set validator
NPOINT_EVAL = (
    cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET
) * cfg.NUM_TIMESTAMPS
residual_validator = ppsci.validate.GeometryValidator(
    equation["NavierStokes"].equations,
    {"continuity": 0, "momentum_x": 0, "momentum_y": 0},
    geom["time_rect_eval"],
    {
        "dataset": "NamedArrayDataset",
        "total_size": NPOINT_EVAL,
        "batch_size": cfg.EVAL.batch_size,
        "sampler": {"name": "BatchSampler"},
    },
    ppsci.loss.MSELoss("mean"),
    metric={"MSE": ppsci.metric.MSE()},
    name="Residual",
)
validator = {residual_validator.name: residual_validator}

The equation setting is the same as the setting of Constraint Construction, indicating how to calculate the target variables to be evaluated;

Here we set the label value to 0 for the three target variables momentum_x, continuity, momentum_y;

The computational domain is the same as the setting of Constraint Construction, indicating evaluation on the specified computational domain;

The sampling point configuration needs to specify the total number of evaluation points total_size. Here we set it to 9662 * 50 (9420 points in the fluid domain + 161 fluid domain inflow boundary points + 81 fluid domain outflow boundary points, a total of 50 evaluation moments);

For evaluation metric metric, select ppsci.metric.MSE;

Other configurations are similar to the settings of Constraint Construction.

3.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 two-dimensional point set in an area. The coordinates of each moment \(t\) are \((x^t_i, y^t_i)\), and the corresponding value is \((u^t_i, v^t_i, p^t_i)\). Therefore, we only need to save the evaluated output data as 50 vtu format files according to time, and finally open them with visualization software to view. The code is as follows:

191
192
193
194
195
196
197
198
199
200
201
202
203
204
# set visualizer(optional)
vis_points = geom["time_rect_eval"].sample_interior(
    (cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET)
    * cfg.NUM_TIMESTAMPS,
    evenly=True,
)
visualizer = {
    "visualize_u_v_p": ppsci.visualize.VisualizerVtu(
        vis_points,
        {"u": lambda d: d["u"], "v": lambda d: d["v"], "p": lambda d: d["p"]},
        num_timestamps=cfg.NUM_TIMESTAMPS,
        prefix="result_u_v_p",
    )
}

3.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.

206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
# initialize solver
solver = ppsci.solver.Solver(
    model,
    constraint,
    cfg.output_dir,
    optimizer,
    None,
    cfg.TRAIN.epochs,
    cfg.TRAIN.iters_per_epoch,
    eval_during_train=cfg.TRAIN.eval_during_train,
    eval_freq=cfg.TRAIN.eval_freq,
    equation=equation,
    geom=geom,
    validator=validator,
    visualizer=visualizer,
    checkpoint_path=cfg.TRAIN.checkpoint_path,
)
# train model
solver.train()
# evaluate after finished training
solver.eval()
# visualize prediction after finished training
solver.visualize()

4. Complete Code

cylinder2d_unsteady_Re100.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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
# 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.

from os import path as osp

import hydra
import numpy as np
from omegaconf import DictConfig

import ppsci
from ppsci.utils import logger
from ppsci.utils import reader


def train(cfg: DictConfig):
    # set random seed for reproducibility
    ppsci.utils.misc.set_random_seed(cfg.seed)

    # initialize logger
    logger.init_logger("ppsci", osp.join(cfg.output_dir, "train.log"), "info")

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

    # set equation
    equation = {
        "NavierStokes": ppsci.equation.NavierStokes(cfg.VISCOSITY, cfg.DENSITY, 2, True)
    }

    # set timestamps
    train_timestamps = np.linspace(
        cfg.TIME_START, cfg.TIME_END, cfg.NUM_TIMESTAMPS, endpoint=True
    ).astype("float32")
    train_timestamps = np.random.choice(train_timestamps, cfg.TRAIN_NUM_TIMESTAMPS)
    train_timestamps.sort()
    t0 = np.array([cfg.TIME_START], dtype="float32")

    val_timestamps = np.linspace(
        cfg.TIME_START, cfg.TIME_END, cfg.NUM_TIMESTAMPS, endpoint=True
    ).astype("float32")

    logger.message(f"train_timestamps: {train_timestamps.tolist()}")
    logger.message(f"val_timestamps: {val_timestamps.tolist()}")

    # set time-geometry
    geom = {
        "time_rect": ppsci.geometry.TimeXGeometry(
            ppsci.geometry.TimeDomain(
                cfg.TIME_START,
                cfg.TIME_END,
                timestamps=np.concatenate((t0, train_timestamps), axis=0),
            ),
            ppsci.geometry.PointCloud(
                reader.load_csv_file(
                    cfg.DOMAIN_TRAIN_PATH,
                    ("x", "y"),
                    alias_dict={"x": "Points:0", "y": "Points:1"},
                ),
                ("x", "y"),
            ),
        ),
        "time_rect_eval": ppsci.geometry.PointCloud(
            reader.load_csv_file(
                cfg.DOMAIN_EVAL_PATH,
                ("t", "x", "y"),
            ),
            ("t", "x", "y"),
        ),
    }

    # pde/bc/sup constraint use t1~tn, initial constraint use t0
    NTIME_PDE = len(train_timestamps)
    ALIAS_DICT = {"x": "Points:0", "y": "Points:1", "u": "U:0", "v": "U:1"}

    # set constraint
    pde_constraint = ppsci.constraint.InteriorConstraint(
        equation["NavierStokes"].equations,
        {"continuity": 0, "momentum_x": 0, "momentum_y": 0},
        geom["time_rect"],
        {
            "dataset": "IterableNamedArrayDataset",
            "batch_size": cfg.NPOINT_PDE * NTIME_PDE,
            "iters_per_epoch": cfg.TRAIN.iters_per_epoch,
        },
        ppsci.loss.MSELoss("mean"),
        name="EQ",
    )
    bc_inlet_cylinder = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "IterableCSVDataset",
                "file_path": cfg.DOMAIN_INLET_CYLINDER_PATH,
                "input_keys": ("x", "y"),
                "label_keys": ("u", "v"),
                "alias_dict": ALIAS_DICT,
                "weight_dict": {"u": 10, "v": 10},
                "timestamps": train_timestamps,
            },
        },
        ppsci.loss.MSELoss("mean"),
        name="BC_inlet_cylinder",
    )
    bc_outlet = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "IterableCSVDataset",
                "file_path": cfg.DOMAIN_OUTLET_PATH,
                "input_keys": ("x", "y"),
                "label_keys": ("p",),
                "alias_dict": ALIAS_DICT,
                "timestamps": train_timestamps,
            },
        },
        ppsci.loss.MSELoss("mean"),
        name="BC_outlet",
    )
    ic = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "IterableCSVDataset",
                "file_path": cfg.IC0_1_PATH,
                "input_keys": ("x", "y"),
                "label_keys": ("u", "v", "p"),
                "alias_dict": ALIAS_DICT,
                "weight_dict": {"u": 10, "v": 10, "p": 10},
                "timestamps": t0,
            },
        },
        ppsci.loss.MSELoss("mean"),
        name="IC",
    )
    sup_constraint = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "IterableCSVDataset",
                "file_path": cfg.PROBE1_50_PATH,
                "input_keys": ("t", "x", "y"),
                "label_keys": ("u", "v"),
                "alias_dict": ALIAS_DICT,
                "weight_dict": {"u": 10, "v": 10},
                "timestamps": train_timestamps,
            },
        },
        ppsci.loss.MSELoss("mean"),
        name="Sup",
    )

    # wrap constraints together
    constraint = {
        pde_constraint.name: pde_constraint,
        bc_inlet_cylinder.name: bc_inlet_cylinder,
        bc_outlet.name: bc_outlet,
        ic.name: ic,
        sup_constraint.name: sup_constraint,
    }

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

    # set validator
    NPOINT_EVAL = (
        cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET
    ) * cfg.NUM_TIMESTAMPS
    residual_validator = ppsci.validate.GeometryValidator(
        equation["NavierStokes"].equations,
        {"continuity": 0, "momentum_x": 0, "momentum_y": 0},
        geom["time_rect_eval"],
        {
            "dataset": "NamedArrayDataset",
            "total_size": NPOINT_EVAL,
            "batch_size": cfg.EVAL.batch_size,
            "sampler": {"name": "BatchSampler"},
        },
        ppsci.loss.MSELoss("mean"),
        metric={"MSE": ppsci.metric.MSE()},
        name="Residual",
    )
    validator = {residual_validator.name: residual_validator}

    # set visualizer(optional)
    vis_points = geom["time_rect_eval"].sample_interior(
        (cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET)
        * cfg.NUM_TIMESTAMPS,
        evenly=True,
    )
    visualizer = {
        "visualize_u_v_p": ppsci.visualize.VisualizerVtu(
            vis_points,
            {"u": lambda d: d["u"], "v": lambda d: d["v"], "p": lambda d: d["p"]},
            num_timestamps=cfg.NUM_TIMESTAMPS,
            prefix="result_u_v_p",
        )
    }

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


def evaluate(cfg: DictConfig):
    # set random seed for reproducibility
    ppsci.utils.misc.set_random_seed(cfg.seed)

    # initialize logger
    logger.init_logger("ppsci", osp.join(cfg.output_dir, "eval.log"), "info")

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

    # set equation
    equation = {
        "NavierStokes": ppsci.equation.NavierStokes(cfg.VISCOSITY, cfg.DENSITY, 2, True)
    }

    # set timestamps
    val_timestamps = np.linspace(
        cfg.TIME_START, cfg.TIME_END, cfg.NUM_TIMESTAMPS, endpoint=True
    ).astype("float32")

    logger.message(f"val_timestamps: {val_timestamps.tolist()}")

    # set time-geometry
    geom = {
        "time_rect_eval": ppsci.geometry.PointCloud(
            reader.load_csv_file(
                cfg.DOMAIN_EVAL_PATH,
                ("t", "x", "y"),
            ),
            ("t", "x", "y"),
        ),
    }

    # set validator
    NPOINT_EVAL = (
        cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET
    ) * cfg.NUM_TIMESTAMPS
    residual_validator = ppsci.validate.GeometryValidator(
        equation["NavierStokes"].equations,
        {"continuity": 0, "momentum_x": 0, "momentum_y": 0},
        geom["time_rect_eval"],
        {
            "dataset": "NamedArrayDataset",
            "total_size": NPOINT_EVAL,
            "batch_size": cfg.EVAL.batch_size,
            "sampler": {"name": "BatchSampler"},
        },
        ppsci.loss.MSELoss("mean"),
        metric={"MSE": ppsci.metric.MSE()},
        name="Residual",
    )
    validator = {residual_validator.name: residual_validator}

    # set visualizer(optional)
    vis_points = geom["time_rect_eval"].sample_interior(
        (cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET)
        * cfg.NUM_TIMESTAMPS,
        evenly=True,
    )
    visualizer = {
        "visualize_u_v_p": ppsci.visualize.VisualizerVtu(
            vis_points,
            {"u": lambda d: d["u"], "v": lambda d: d["v"], "p": lambda d: d["p"]},
            num_timestamps=cfg.NUM_TIMESTAMPS,
            prefix="result_u_v_p",
        )
    }

    # initialize solver
    solver = ppsci.solver.Solver(
        model,
        geom=geom,
        output_dir=cfg.output_dir,
        validator=validator,
        visualizer=visualizer,
        pretrained_model_path=cfg.EVAL.pretrained_model_path,
    )
    # evaluate
    solver.eval()
    # visualize prediction
    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 time-geometry
    geom = {
        "time_rect_eval": ppsci.geometry.PointCloud(
            reader.load_csv_file(
                cfg.DOMAIN_EVAL_PATH,
                ("t", "x", "y"),
            ),
            ("t", "x", "y"),
        ),
    }
    NPOINT_EVAL = (
        cfg.NPOINT_PDE + cfg.NPOINT_INLET_CYLINDER + cfg.NPOINT_OUTLET
    ) * cfg.NUM_TIMESTAMPS
    input_dict = geom["time_rect_eval"].sample_interior(NPOINT_EVAL, evenly=True)
    output_dict = predictor.predict(input_dict, 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())
    }

    ppsci.visualize.save_vtu_from_dict(
        "./cylinder2d_unsteady_Re100_pred.vtu",
        {**input_dict, **output_dict},
        input_dict.keys(),
        cfg.MODEL.output_keys,
        cfg.NUM_TIMESTAMPS,
    )


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

The prediction results are displayed below. The horizontal axis represents the x-direction, and the vertical axis represents the y-direction, with fluid flowing from left to right. The figure illustrates the predicted lateral flow velocity \(u(t,x,y)\) at 50 time steps.

Note

This case is only shown as a demo and has not been fully tuned. Some of the results shown below may differ from OpenFOAM.

u_pred.gif

Model prediction result of lateral flow velocity u