关于三类边界条件 #1083
Comments
|
你好，PaddleScience里并没有直接实现某一种类型的边界条件，而是希望用户基于自己的公式手动编写，从而支持任意类型的边界条件。 x = sympy.symbols('x')
u = sympy.Function('u')(x)
def du_dx_pred(data):
# 注意这里的data['h']和data['u0']是你需要在dataloader中准备好的输入数据，data['u']是模型本身的输出，会自动添加到data中
return data['h'] * (data['u'] - data['u0'])
bc_inlet = ppsci.constraint.BoundaryConstraint(
{"du/dx": u.diff(x)}, # 等式左侧公式
{"du/dx": du_dx_pred}, # 等式右侧计算公式
your_geometry_object,
{**train_dataloader_cfg, "batch_size": cfg.TRAIN.batch_size.bc_inlet},
ppsci.loss.MSELoss("sum"),
name="inlet",
)我们有很多的带有复杂边界条件的案例可以参考，如：https://paddlescience-docs.readthedocs.io/zh-cn/latest/zh/examples/aneurysm/#342 |
|
如果一开始方程的定义是用的ppsci里带的方程，并没有自己去使用sympy定义输入输出变量，比如equation = {"laplace": ppsci.equation.Laplace(dim=2)}}，边界条件有导数，比如du/dx, du/dy，那怎么做呢？ 是再定义一遍 还是能够下面这样？ |
测试了一下，如果这样就可以。不是很理解原理 |
|
@pecanjk ppsci中的XXXConstraint有两个关键参数: PaddleScience/ppsci/utils/expression.py Lines 164 to 166 in d58a744 如果是sympy表达式，那么就会由 PaddleScience/ppsci/solver/solver.py Lines 507 to 546 in d58a744 然后再在expression.py中被同样的方式调用 你写的那份能跑通的代码其实是不正确的，因为你写的u是一个sympy变量对象，而不是一个sympy函数对象，所以u.diff(x)其实结果是0，而不是表达式，所以你的代码仅仅能跑通，而并不是正确的。 需要自定义方程，请参考我们的文档：https://paddlescience-docs.readthedocs.io/zh-cn/latest/zh/development/#241，【2.4 构建方程】章节 |
|
方程自定义后，还是没法拿到输出变量的一阶导，除非自定义方程中直接把一阶导数也用 我看了nvidia modulus的案例的modulus-sym/examples/wave_equation/wave_1d.py，对于wave equation, du/dt(0)一阶导的初值边界条件会经常用，它提供了'u__t'这样的计算 |
输出变量的一阶导数是如何计算的，需要用户指定，modulus里用的是双下划线 |
|
或者你可以提供你自己的写的方程代码，我可以帮你看一下 |
|
您看一下，wave equation ##main.py
import hydra
import numpy as np
from omegaconf import DictConfig
import paddle
import ppsci
from wave_equation import Wave
c=1
# compute ground truth function
def u_solution_func(data):
"""compute ground truth for u as label data"""
t, x = data["t"], data["x"]
return np.sin(np.pi * x) * np.cos(c * np.pi * t) + np.sin(2 * np.pi * x) * np.cos(2 * c * np.pi * t)
def train(cfg: DictConfig):
# set random seed for reproducibility
ppsci.utils.misc.set_random_seed(cfg.seed)
# set equation
equation = {"wave": Wave(dim=1, time=True, c=c)}
# set geometry domain
t_domain = ppsci.geometry.TimeDomain(**cfg.TIME_COORD) # 时间域
geom_domain = ppsci.geometry.Interval(**cfg.GEOM_COORD)
domain = ppsci.geometry.TimeXGeometry(t_domain, geom_domain)
# set model
model = ppsci.arch.MLP(**cfg.MODEL)
train_dataloader_cfg = {
"dataset": "IterableNamedArrayDataset",
"iters_per_epoch": cfg.TRAIN.iters_per_epoch,
}
# set constraint
pde_constraint = ppsci.constraint.InteriorConstraint(
equation["wave"].equations,
{"wave": 0},
domain,
{**train_dataloader_cfg, "batch_size": cfg.NPOINT_TOTAL},
ppsci.loss.CausalMSELoss(20,"mean",2),
evenly=True,
name="EQ",
)
bc = ppsci.constraint.BoundaryConstraint(
{"u": lambda out: out["u"]},
{"u": 0},
domain,
{**train_dataloader_cfg, "batch_size": cfg.NPOINT_BC},
ppsci.loss.MSELoss("mean"),
evenly=True,
name="BC",
)
ic = ppsci.constraint.InitialConstraint(
{"u0": lambda out: out["u"], "u_t0": u.diff(t)}, ##这里需要一阶导的计算
{"u0": lambda _in: np.sin(np.pi * _in["x"]) + np.sin(2* np.pi * _in["x"]), "u_t0": 0},
domain,
{**train_dataloader_cfg, "batch_size": cfg.NPOINT_IC},
ppsci.loss.MSELoss("mean"),
evenly=True,
name="IC",
)
# wrap constraints together
constraint = {
pde_constraint.name: pde_constraint,
bc.name: bc,
ic.name: ic,
}
# set optimizer
optimizer = ppsci.optimizer.Adam(learning_rate=cfg.TRAIN.learning_rate)(model)
# set validator
mse_metric = ppsci.validate.GeometryValidator(
{"u": lambda out: out["u"]},
{"u": u_solution_func},
domain,
{
"dataset": "IterableNamedArrayDataset",
"total_size": cfg.NPOINT_TOTAL,
},
ppsci.loss.MSELoss(),
evenly=True,
metric={"MSE": ppsci.metric.MSE()},
with_initial=False,
name="MSE_Metric",
)
validator = {mse_metric.name: mse_metric}
# initialize solver
solver = ppsci.solver.Solver(
model,
constraint,
output_dir=cfg.output_dir,
optimizer=optimizer,
validator=validator,
cfg=cfg
)
# train model
solver.train()
@hydra.main(version_base=None, config_path="./conf", config_name="wave_1d.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()##自定义的wave_equation
from __future__ import annotations
from typing import Optional
from typing import Tuple
from ppsci.equation import PDE
class Wave(PDE):
def __init__(self, dim: int=1, time: bool=True, c = 1.0, detach_keys: Optional[Tuple[str, ...]] = None):
super().__init__()
self.detach_keys = detach_keys
self.dim = dim
self.time = time
t, x, y, z = self.create_symbols("t x y z")
if dim == 1:
invars = (x,)
elif dim == 2:
invars = (x, y)
elif dim == 3:
invars = (x, y, z)
else:
raise ValueError("Dimension should be 1, 2 or 3.")
if time:
invars = (t,) + invars
u = self.create_function("u", invars)
if time:
u_tt = u.diff(t,2)
coord_var = invars[1:]
else:
u_tt = 0
coord_var = invars
poisson = 0
for invar in coord_var:
poisson += u.diff(invar, 2)
self.add_equation("wave", u_tt - c**2 * poisson)
self._apply_detach()defaults:
- ppsci_default
- TRAIN: train_default
- TRAIN/ema: ema_default
- TRAIN/swa: swa_default
- EVAL: eval_default
- INFER: infer_default
- hydra/job/config/override_dirname/exclude_keys: exclude_keys_default
- _self_
hydra:
run:
# dynamic output directory according to running time and override name
dir: outputs_wave_1d/${now:%Y-%m-%d}/${now:%H-%M-%S}/${hydra.job.override_dirname}
job:
name: ${mode} # name of logfile
chdir: false # keep current working directory unchanged
callbacks:
init_callback:
_target_: ppsci.utils.callbacks.InitCallback
sweep:
# output directory for multirun
dir: ${hydra.run.dir}
subdir: ./
# general settings
mode: train # running mode: train/eval
seed: 201
output_dir: ${hydra:run.dir}
log_freq: 20
# set geometry
GEOM_COORD:
l: 0
r: 1
TIME_COORD:
t0: 0.0
t1: 1.0
time_step: 0.05
NPOINT_TOTAL: 20000
NPOINT_BC: 100
NPOINT_IC: 200
# model settings
MODEL:
input_keys: ["t", "x"]
output_keys: ["u"]
num_layers: 3
hidden_size: 60
activation: "tanh"
# training settings
TRAIN:
epochs: 4000
iters_per_epoch: 1
eval_during_train: true
eval_freq: 100
learning_rate: 1.0e-3
pretrained_model_path: null
# evaluation settings
EVAL:
pretrained_model_path: null
# inference settings
INFER:
pretrained_model_path: null |
|
看到你的方程仅定义在Wave中，所以你的sympy函数、变量也只有在Wave这个PDE中生效，所以可以在main.py中再次定义函数u、变量x、t，然后使用；或者用ppsci.autodiff下的两个接口也可以 ## main.py
### 第一种方式
t, x = sympy.symbols('t x')
u = sympy.Function('u')(t, x)
ic = ppsci.constraint.InitialConstraint(
{"u0": lambda out: out["u"], "u_t0": u.diff(t)},
{"u0": lambda _in: np.sin(np.pi * _in["x"]) + np.sin(2* np.pi * _in["x"]), "u_t0": 0},
domain,
{**train_dataloader_cfg, "batch_size": cfg.NPOINT_IC},
ppsci.loss.MSELoss("mean"),
evenly=True,
name="IC",
)
### 第二种方式
from ppsci.autodiff import jacobian, hessian
ic = ppsci.constraint.InitialConstraint(
{"u0": lambda out: out["u"], "u_t0": lambda out: jacobian(out["u"], out["t"])},
{"u0": lambda _in: np.sin(np.pi * _in["x"]) + np.sin(2* np.pi * _in["x"]), "u_t0": 0},
domain,
{**train_dataloader_cfg, "batch_size": cfg.NPOINT_IC},
ppsci.loss.MSELoss("mean"),
evenly=True,
name="IC",
) |
|
感谢，能否加入官方案例中去，目前还没有wave equation的支持 |
加入官方案例的话需要开发者提交PR，包含代码、文档、训练权重，可以参考：https://paddlescience-docs.readthedocs.io/zh-cn/latest/zh/development/#3。如果你有时间和兴趣的话，可以把你的这个案例完善一下，然后提交PR到这个仓库里，我会帮忙review |
|
非常感谢讲解 @HydrogenSulfate ，现在我有点了解应该怎么写边界条件这部分了。如果我的约束条件中的u也是网络输出，而我希望约束du/dx，是不是能像下面这么写： 其中h和u0是定义好的常量。 |
|
你理解的是正确的，但是需要注意一下sympy、Tensor两种不同类型下的微分写法，out是
改为： Tensor和sympy两种写法虽然可以同时共存使用，但这是两套体系，是不能相互混用对方的函数 然后下面这行没问题
|
|
@HydrogenSulfate 像上面那样写会报错： 其中equation = {"heat_pde": ppsci.equation.HeatTransfer(alpha=45.0, dim=3)}这个是我自己的方程 |
|
@ZBDat 在一个三维的立方体geom表面随机撒点，筛选出足够多的y=1其实是比较困难的，所以建议你的bc可以再重新构建一个y=1附近的三维的薄片（仍然可以使用Cuboid，把对角线y的上下界设置的非常接近1即可），然后把这个三维薄片作为bc的几何即可 |
请提出你的问题 Please ask your question
请问边界条件类是否支持neumann边界条件？第三类边界条件呢？比如du/dx = h (u - u0)
The text was updated successfully, but these errors were encountered: