{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 基于MindFlow定义符号化偏微分方程\n",
    "\n",
    "[![下载Notebook](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.0.0-alpha/resource/_static/logo_notebook.png)](https://obs.dualstack.cn-north-4.myhuaweicloud.com/mindspore-website/notebook/r2.0.0-alpha/mindflow/zh_cn/physics_driven/mindspore_sympy_pde_definition.ipynb) [![下载样例代码](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.0.0-alpha/resource/_static/logo_download_code.png)](https://obs.dualstack.cn-north-4.myhuaweicloud.com/mindspore-website/notebook/r2.0.0-alpha/mindflow/zh_cn/physics_driven/mindspore_sympy_pde_definition.py) [![查看源文件](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/r2.0.0-alpha/resource/_static/logo_source.png)](https://gitee.com/mindspore/docs/blob/r2.0.0-alpha/docs/mindflow/docs/source_zh_cn/physics_driven/sympy_pde_definition.ipynb)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "偏微分方程（PDE）在工程应用中发挥着重要作用，大多数控制系统和人造复杂系统的物理学都是有PDE描述的。本文介绍如何基于MindFlow使用sympy定义PDE，构建基于物理神经网络（Physics-Informed Neural Networks，PINNs）的问题求解模型。利用 *mindflow.pde.PDEWithLoss* ，可以很容易以符号形式描绘偏微分方程，并计算所有方程的损失。这使得方程的描述简洁易懂，易于扩展。用户可继承 *mindflow.pde.PDEWithLoss* 实现自定义偏微分方程。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 偏微分方程样例：Navier-Stokes equation\n",
    "\n",
    "Navier-Stokes方程，被称为 `N-S` 方程，是流体力学领域的经典偏微分方程。在粘性不可压缩的情况下，无量纲 `N-S` 方程具有以下形式：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y} = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u} {\\partial t} + u \\frac{\\partial u}{\\partial x} + v \\frac{\\partial u}{\\partial y} + \\frac{\\partial p}{\\partial x} - \\frac{1} {Re} (\\frac{\\partial^2u}{\\partial x^2} + \\frac{\\partial^2u}{\\partial y^2}) = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial v} {\\partial t} + u \\frac{\\partial v}{\\partial x} + v \\frac{\\partial v}{\\partial y} + \\frac{\\partial p}{\\partial y} - \\frac{1} {Re} (\\frac{\\partial^2v}{\\partial x^2} + \\frac{\\partial^2v}{\\partial y^2}) = 0\n",
    "$$\n",
    "\n",
    "其中 `Re` 表示雷诺数。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 符号声明"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[x, y, t]\n",
      "[u(x, y, t), v(x, y, t), p(x, y, t)]\n"
     ]
    }
   ],
   "source": [
    "from sympy import symbols, Function\n",
    "\n",
    "x, y, t = symbols('x y t')\n",
    "u = Function('u')(x, y, t)\n",
    "v = Function('v')(x, y, t)\n",
    "p = Function('p')(x, y, t)\n",
    "\n",
    "# independent variables\n",
    "in_vars = [x, y, t]\n",
    "print(in_vars)\n",
    "\n",
    "# dependent variables\n",
    "out_vars = [u, v, p]\n",
    "print(out_vars)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 符号形式的偏微分方程\n",
    "\n",
    "使用上述符号声明定义偏微分方程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sympy import diff"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 控制方程"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u(x, y, t)*Derivative(u(x, y, t), x) + v(x, y, t)*Derivative(u(x, y, t), y) + Derivative(p(x, y, t), x) + Derivative(u(x, y, t), t) - 0.00999999977648258*Derivative(u(x, y, t), (x, 2)) - 0.00999999977648258*Derivative(u(x, y, t), (y, 2))\n",
      "u(x, y, t)*Derivative(v(x, y, t), x) + v(x, y, t)*Derivative(v(x, y, t), y) + Derivative(p(x, y, t), y) + Derivative(v(x, y, t), t) - 0.00999999977648258*Derivative(v(x, y, t), (x, 2)) - 0.00999999977648258*Derivative(v(x, y, t), (y, 2))\n",
      "Derivative(u(x, y, t), x) + Derivative(v(x, y, t), y)\n"
     ]
    }
   ],
   "source": [
    "# Consider Reynolds number is 100\n",
    "re = 100\n",
    "number = np.float32(1.0 / re)\n",
    "\n",
    "# X Momemtum\n",
    "momentum_x = u.diff(t) + u * u.diff(x) + v * u.diff(y) + p.diff(x) - number*(diff(u, (x, 2)) + diff(u, (y, 2)))\n",
    "print(momentum_x)\n",
    "\n",
    "# Y Momemtum\n",
    "momentum_y = v.diff(t) + u * v.diff(x) + v * v.diff(y) + p.diff(y) - number*(diff(v, (x, 2)) + diff(v, (y, 2)))\n",
    "print(momentum_y)\n",
    "\n",
    "# continuty\n",
    "continuty = u.diff(x) + v.diff(y)\n",
    "print(continuty)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 边界条件\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u(x, y, t)\n",
      "v(x, y, t)\n"
     ]
    }
   ],
   "source": [
    "bc_u = u\n",
    "print(bc_u)\n",
    "\n",
    "bc_v = v\n",
    "print(bc_v)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 初始条件"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u(x, y, t)\n",
      "v(x, y, t)\n",
      "p(x, y, t)\n"
     ]
    }
   ],
   "source": [
    "ic_u = u\n",
    "print(bc_u)\n",
    "\n",
    "ic_v = v\n",
    "print(ic_v)\n",
    "\n",
    "ic_p = p\n",
    "print(ic_p)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题建模样例\n",
    "\n",
    "下面的 `CylinderFlow` 定义了二维非定常圆柱绕流问题。具体来说，它包括上面定义的三个部分：控制方程、初始条件和边界条件。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`NavierStokes` 基类，定义了输入和输出变量，以及控制方程。\n",
    "\n",
    "用户可以通过重写 `NavierStokes` 类定义其他控制方程。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mindspore import nn\n",
    "from mindflow.pde import PDEWithLoss\n",
    "\n",
    "\n",
    "class NavierStokes(PDEWithLoss):\n",
    "    def __init__(self, model, re=100, loss_fn=nn.MSELoss()):\n",
    "        self.number = np.float32(1.0 / re)\n",
    "        self.x, self.y, self.t = symbols('x y t')\n",
    "        self.u = Function('u')(self.x, self.y, self.t)\n",
    "        self.v = Function('v')(self.x, self.y, self.t)\n",
    "        self.p = Function('p')(self.x, self.y, self.t)\n",
    "        self.in_vars = [self.x, self.y, self.t]\n",
    "        self.out_vars = [self.u, self.v, self.p]\n",
    "        super(NavierStokes, self).__init__(model, self.in_vars, self.out_vars)\n",
    "        self.loss_fn = loss_fn\n",
    "\n",
    "    def pde(self):\n",
    "        momentum_x = self.u.diff(self.t) + self.u * self.u.diff(self.x) + self.v * self.u.diff(self.y) +\\\n",
    "            self.p.diff(self.x) - self.number * (diff(self.u, (self.x, 2)) + diff(self.u, (self.y, 2)))\n",
    "        momentum_y = self.v.diff(self.t) + self.u * self.v.diff(self.x) + self.v * self.v.diff(self.y) +\\\n",
    "            self.p.diff(self.y) - self.number * (diff(self.v, (self.x, 2)) + diff(self.v, (self.y, 2)))\n",
    "        continuty = self.u.diff(self.x) + self.v.diff(self.y)\n",
    "\n",
    "        equations = {\"momentum_x\": momentum_x, \"momentum_y\": momentum_y, \"continuty\": continuty}\n",
    "        return equations"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面，我们使用 `NavierStokes` 基类定义初始条件和边界条件，以及损失函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mindspore import nn, ops, Tensor\n",
    "from mindspore import dtype as mstype\n",
    "from mindflow.pde import PDEWithLoss, NavierStokes, sympy_to_mindspore\n",
    "\n",
    "\n",
    "class CylinderFlow(NavierStokes):\n",
    "    def __init__(self, model, re=100, loss_fn=nn.MSELoss()):\n",
    "        super(CylinderFlow, self).__init__(model, re=re, loss_fn=loss_fn)\n",
    "        self.ic_nodes = sympy_to_mindspore(self.ic(), self.in_vars, self.out_vars)\n",
    "        self.bc_nodes = sympy_to_mindspore(self.bc(), self.in_vars, self.out_vars)\n",
    "\n",
    "    def bc(self):\n",
    "        bc_u = self.u\n",
    "        bc_v = self.v\n",
    "        equations = {\"bc_u\": bc_u, \"bc_v\": bc_v}\n",
    "        return equations\n",
    "\n",
    "    def ic(self):\n",
    "        ic_u = self.u\n",
    "        ic_v = self.v\n",
    "        ic_p = self.p\n",
    "        equations = {\"ic_u\": ic_u, \"ic_v\": ic_v, \"ic_p\": ic_p}\n",
    "        return equations\n",
    "\n",
    "    def get_loss(self, pde_data, bc_data, bc_label, ic_data, ic_label):\n",
    "        pde_res = self.parse_node(self.pde_nodes, inputs=pde_data)\n",
    "        pde_residual = ops.Concat(1)(pde_res)\n",
    "        pde_loss = self.loss_fn(pde_residual, Tensor(np.array([0.0]).astype(np.float32), mstype.float32))\n",
    "\n",
    "        ic_res = self.parse_node(self.ic_nodes, inputs=ic_data)\n",
    "        ic_residual = ops.Concat(1)(ic_res)\n",
    "        ic_loss = self.loss_fn(ic_residual, ic_label)\n",
    "\n",
    "        bc_res = self.parse_node(self.bc_nodes, inputs=bc_data)\n",
    "        bc_residual = ops.Concat(1)(bc_res)\n",
    "        bc_loss = self.loss_fn(bc_residual, bc_label)\n",
    "\n",
    "        return pde_loss + ic_loss + bc_loss\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "详细案例参见 [基于PINNs关于圆柱绕流的Navier-Stokes equation求解](https://www.mindspore.cn/mindflow/docs/zh-CN/r0.1.0-alpha/physics_driven/navier_stokes2D.html)。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neumann边界条件定义样例"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在数学中，Neumann边界条件也被称为常微分方程或偏微分方程的“第二类边界条件”。Neumann边界条件指定微分方程解边界处的微分。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下述 `Poisson2D` 问题定义了Dirichlet边界条件(`bc_outer`)和Neumann边界条件(`bc_inner`)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "from mindflow.pde import Poisson\n",
    "\n",
    "class Poisson2D(Poisson):\n",
    "    def __init__(self, model, loss_fn=nn.MSELoss()):\n",
    "        super(Poisson2D, self).__init__(model, loss_fn=loss_fn)\n",
    "        self.bc_outer_nodes = sympy_to_mindspore(self.bc_outer(), self.in_vars, self.out_vars)\n",
    "        self.bc_inner_nodes = sympy_to_mindspore(self.bc_inner(), self.in_vars, self.out_vars)\n",
    "\n",
    "    def bc_outer(self):\n",
    "        bc_outer_eq = self.u\n",
    "        equations = {\"bc_outer\": bc_outer_eq}\n",
    "        return equations\n",
    "\n",
    "    def bc_inner(self):\n",
    "        bc_inner_eq = sympy.Derivative(self.u, self.normal) - 0.5\n",
    "        equations = {\"bc_inner\": bc_inner_eq}\n",
    "        return equations\n",
    "\n",
    "    def get_loss(self, pde_data, bc_outer_data, bc_inner_data, bc_inner_normal):\n",
    "        pde_res = self.parse_node(self.pde_nodes, inputs=pde_data)\n",
    "        pde_loss = self.loss_fn(pde_res[0], Tensor(np.array([0.0]), mstype.float32))\n",
    "\n",
    "        bc_inner_res = self.parse_node(self.bc_inner_nodes, inputs=bc_inner_data, norm=bc_inner_normal)\n",
    "        bc_inner_loss = self.loss_fn(bc_inner_res[0], Tensor(np.array([0.0]), mstype.float32))\n",
    "\n",
    "        bc_outer_res = self.parse_node(self.bc_outer_nodes, inputs=bc_outer_data)\n",
    "        bc_outer_loss = self.loss_fn(bc_outer_res[0], Tensor(np.array([0.0]), mstype.float32))\n",
    "\n",
    "        return pde_loss + bc_inner_loss + bc_outer_loss"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "详细案例参见 [基于PINNs关于作用于圆环的Poisson equation求解](https://www.mindspore.cn/mindflow/docs/zh-CN/r0.1.0-alpha/physics_driven/poisson_ring.html)。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.0 ('py39')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.0"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "57ace93c29d9374277a79956c3f1b916d7d9a05468d906842f9921d0d494a29f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}