cleanup of forw ns2

This commit is contained in:
NT
2021-02-18 16:00:22 +08:00
parent 56dc78c32c
commit 8e0404b259
2 changed files with 86 additions and 31 deletions

View File

@@ -8,6 +8,11 @@
"source": [
"# Differentiable Fluid Simulations\n",
"\n",
"\n",
"**TODO , NEXT update diffphys example, then check burgers examples**\n",
"\n",
"\n",
"\n",
"Next, we'll target a more complex example with the Navier-Stokes equations as model. We'll target a 2D case with velocity $\\mathbf{u}$, no explicit viscosity term, and a marker density $d$ that drives a simple Boussinesq buoyancy term $\\eta d$ adding a force along the y dimension:\n",
"\n",
"$\\begin{aligned}\n",

View File

@@ -8,8 +8,6 @@
"source": [
"# Navier-Stokes Forward Simulation\n",
"\n",
"**TODO , double check here - sync with latest local-colab-notebook version, then update diffphys example**\n",
"\n",
"Now let's target a somewhat more complex example: a fluid simulations based on the Navier-Stokes equations. This is still very simple with Φ<sub>Flow</sub>, as differentiable operators for all steps exist. The Navier-Stokes equations (in their incompressible form) introduce an additional pressure field $p$, and a constraint for conservation of mass:\n",
"\n",
"$\\begin{align}\n",
@@ -69,14 +67,14 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 2,
"metadata": {
"id": "WrA3IXDxv31P"
},
"outputs": [],
"source": [
"DOMAIN = Domain(x=40, y=32, boundaries=CLOSED, bounds=Box[0:100, 0:80])\n",
"INFLOW = domain.grid(Sphere(center=(30, 15), radius=10)) * 0.2\n",
"INFLOW = DOMAIN.grid(Sphere(center=(30, 15), radius=10)) * 0.2\n",
"\n",
"DT = 1.5\n",
"NU = 0.01"
@@ -98,7 +96,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
@@ -118,7 +116,7 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
@@ -138,10 +136,10 @@
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7f7f88ebd4c0>"
"<matplotlib.image.AxesImage at 0x7fb2c07a7c40>"
]
},
"execution_count": 12,
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
@@ -159,15 +157,15 @@
}
],
"source": [
"def step(velocity, smoke, pressure, buoyancy_factor=1.0):\n",
" smoke = advect.semi_lagrangian(smoke, velocity, dt) + inflow\n",
"def step(velocity, smoke, pressure, dt=1.0, buoyancy_factor=1.0):\n",
" smoke = advect.semi_lagrangian(smoke, velocity, dt) + INFLOW\n",
" buoyancy_force = smoke * (0, buoyancy_factor) >> velocity # resamples smoke to velocity sample points\n",
" velocity = advect.semi_lagrangian(velocity, velocity, 1) + dt * buoyancy_force\n",
" velocity = diffuse.explicit(velocity, nu, dt)\n",
" velocity, pressure, iterations, divergence = fluid.make_incompressible(velocity, domain, pressure_guess=pressure)\n",
" velocity = diffuse.explicit(velocity, NU, dt)\n",
" velocity, pressure, iterations, divergence = fluid.make_incompressible(velocity, DOMAIN, pressure_guess=pressure)\n",
" return velocity, smoke, pressure\n",
"\n",
"velocity, smoke, pressure = step(velocity, smoke, None)\n",
"velocity, smoke, pressure = step(velocity, smoke, None, dt=DT)\n",
"\n",
"print(\"Max. velocity and mean density: \" + format( [ math.max(velocity.values) , math.mean(smoke.values) ] ))\n",
"\n",
@@ -184,36 +182,88 @@
"\n",
"The pressure projection step in `make_incompressible` is typically the computationally most expensive step in the sequence above. It solves a Poisson equation for the boundary conditions of the domain, and updates the velocity field with the gradient of the computed pressure.\n",
"\n",
"Just for testing, we've also printed the mean value of the velocities, and the max density after the update. As you can see in the resulting image, we have a first round region of smoke, with a slight upwards motion (which does not show here yet). "
"Just for testing, we've also printed the mean value of the velocities, and the max density after the update. As you can see in the resulting image, we have a first round region of smoke, with a slight upwards motion (which does not show here yet). \n",
"\n",
"## Datatypes and Dimensions\n",
"\n",
"The created grids are instances of the class `Grid`.\n",
"Like tensors, grids also have the `shape` attribute which lists all batch, spatial and channel dimensions.\n",
"[Shapes in Φ<sub>Flow</sub>](https://tum-pbs.github.io/PhiFlow/Math.html#shapes) store not only the sizes of the dimensions but also their names and types."
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Smoke: (x=40, y=32)\n",
"Velocity: (x=40, y=32, vector=2)\n",
"Inflow: (x=40, y=32)\n",
"Inflow, spatial only: (x=40, y=32)\n"
]
}
],
"source": [
"# ? print(\"Fluid state: \" + format(fluid.state))\n",
"# ? print(\"\\nVelocity content:\")\n",
"# ? [print(grid.data.shape) for grid in fluid.velocity.unstack()];"
"print(f\"Smoke: {smoke.shape}\")\n",
"print(f\"Velocity: {velocity.shape}\")\n",
"print(f\"Inflow: {INFLOW.shape}\")\n",
"print(f\"Inflow, spatial only: {INFLOW.shape.spatial}\")\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TODO\n",
"\n",
"Note that we actually created two variables, one for each velocity component. If you're interested in how this works, have a look at the [Struct documentation](https://github.com/tum-pbs/PhiFlow/blob/master/documentation/Structs.ipynb).\n",
"\n",
"If you look closely at the output from the last `print` command, you'll notice that the shapes of the variables differ. This is because the velocity is sampled in [staggered form](https://github.com/tum-pbs/PhiFlow/blob/master/documentation/Staggered_Grids.md).\n",
"The grid values can be accessed using the `values` property."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(x=40, y=32) float32 0.0 < ... < 0.2\n",
"(x=(41, 40) along vector, y=(32, 33) along vector, vector=2) float32 -0.10927753 < ... < 0.135354\n",
"(x=40, y=32) float32 0.0 < ... < 0.2\n"
]
}
],
"source": [
"print(smoke.values)\n",
"print(velocity.values)\n",
"print(INFLOW.values)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Grids have many more properties which are documented [here](https://tum-pbs.github.io/PhiFlow/phi/field/#phi.field.Grid).\n",
"Also note that the staggered grid has a [non-uniform shape](https://tum-pbs.github.io/PhiFlow/Math.html#non-uniform-tensors) because the number of faces is not equal to the number of cells. The `INFLOW` grid naturally has the same dimensions as the `smoke` grid.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Time Evolution\n",
"\n",
"With this setup, we can easily advance the simulation forward in time a bit more by repeatedly calling the `step` function."
]
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 7,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
@@ -241,7 +291,7 @@
],
"source": [
"for time_step in range(10):\n",
" velocity, smoke, pressure = step(velocity, smoke, pressure)\n",
" velocity, smoke, pressure = step(velocity, smoke, pressure, dt=DT)\n",
" print('Computed frame {}, max velocity {}'.format(time_step , np.asarray(math.max(velocity.values)) ))\n"
]
},
@@ -256,7 +306,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 8,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
@@ -269,10 +319,10 @@
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7f7f8390fdf0>"
"<matplotlib.image.AxesImage at 0x7fb2bc1153d0>"
]
},
"execution_count": 6,
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
@@ -304,7 +354,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
@@ -342,7 +392,7 @@
"for time_step in range(20):\n",
" if time_step<3 or time_step%10==0: \n",
" print('Computing time step %d' % time_step)\n",
" velocity, smoke, pressure = step(velocity, smoke, pressure)\n",
" velocity, smoke, pressure = step(velocity, smoke, pressure, dt=DT)\n",
" if time_step%5==0:\n",
" steps.append(smoke.values.numpy('y,x'))\n",
"\n",