From 1ed0e3ab7b4881b48b76751b0224ae31f471731f Mon Sep 17 00:00:00 2001 From: BillMaZengou Date: Thu, 22 Feb 2024 21:56:33 +0800 Subject: [PATCH] add advdiffreac.ipynb for demo --- chapter2/advdiffreac.ipynb | 337 +++++++++++++++++++++++++++++++++++++ 1 file changed, 337 insertions(+) create mode 100644 chapter2/advdiffreac.ipynb diff --git a/chapter2/advdiffreac.ipynb b/chapter2/advdiffreac.ipynb new file mode 100644 index 00000000..562e187e --- /dev/null +++ b/chapter2/advdiffreac.ipynb @@ -0,0 +1,337 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1b721ef5-7000-4a6c-898f-0ce3cd65a2b0", + "metadata": {}, + "source": [ + "# A system of advection-diffusion-reaction equations\n", + "\n", + "$$\n", + "\\frac{\\partial u_1}{\\partial t} + w \\cdot \\nabla u_1 - \\nabla \\cdot (\\epsilon \\nabla u_1) = f_1 - K u_1 u_2\n", + "$$\n", + "$$\n", + "\\frac{\\partial u_2}{\\partial t} + w \\cdot \\nabla u_2 - \\nabla \\cdot (\\epsilon \\nabla u_2) = f_2 - K u_1 u_2\n", + "$$\n", + "$$\n", + "\\frac{\\partial u_3}{\\partial t} + w \\cdot \\nabla u_3 - \\nabla \\cdot (\\epsilon \\nabla u_3) = f_3 + K u_1 u_2 - K u_3\n", + "$$\n", + "\n", + "This system models the chemical reaction between two species A and B in some domain $\\Omega$\n", + "$$\n", + "A + B \\rightarrow C\n", + "$$\n", + "We assume that the reaction is *first-order*, meaning that the reaction rate is proportional to the concentrations $[A]$ and $[B]$ of the two species:\n", + "$$\n", + "\\frac{d}{dt} [C] = K[A][B]\n", + "$$\n", + "\n", + "$$\n", + "u_1 = [A], u_2 = [B], u_3 = [C]\n", + "$$\n", + "\n", + "The chemical reactions take part at each point in the domain $\\Omega$. In addition, we assume that the species $A, B, C$ diffuse throughout the domain with **diffusivity** $\\epsilon$ and are advected with velocity $w$.\n", + "\n", + "To make things visually and physically interesting, we shall let the chemical reaction take place in the velocity field computed from the solution of imcompressible NS Eq. around a cylinder. i.e.)\n", + "$$\n", + "\\rho \\left(\\frac{\\partial w}{\\partial t} + w \\cdot \\nabla w \\right) = \\nabla \\cdot \\sigma(w, p) + f\n", + "$$\n", + "$$\n", + "\\nabla \\cdot w = 0\n", + "$$\n", + "$$\n", + "\\frac{\\partial u_1}{\\partial t} + w \\cdot \\nabla u_1 - \\nabla \\cdot (\\epsilon \\nabla u_1) = f_1 - K u_1 u_2\n", + "$$\n", + "$$\n", + "\\frac{\\partial u_2}{\\partial t} + w \\cdot \\nabla u_2 - \\nabla \\cdot (\\epsilon \\nabla u_2) = f_2 - K u_1 u_2\n", + "$$\n", + "$$\n", + "\\frac{\\partial u_3}{\\partial t} + w \\cdot \\nabla u_3 - \\nabla \\cdot (\\epsilon \\nabla u_3) = f_3 + K u_1 u_2 - K u_3\n", + "$$\n", + "\n", + "We assume that $u_1 = u_2 = u_3 = 0$ at $t = 0$ and inject the species $A, B$ into the system by specifying non-zero source terms $f_1, f_2$ close to the corners at the inflow, and $f_3 = 0$. The result will be that $A, B$ are convected by advection and diffusion throughout the channel, and when they mix the species $C$ will be formed.\n", + "\n", + "Since the system is one-way coupled from the NS subsystem to the advection-diffusion-reaction subsystem, we do not need to recompute the solution to NS Eq., but can just **read back** the previously computed velocity field and **feed into** out equations. (Learn how to read and write solutions)" + ] + }, + { + "cell_type": "markdown", + "id": "29a53fa3-23e8-40bf-a9e9-93bcc5510fbb", + "metadata": {}, + "source": [ + "## Variational Formulation\n", + "We obtain the variational formulation of our system by multiplying each equation by a test function, integrating the second-order terms $-\\nabla \\cdot (\\epsilon \\nabla u_i)$ by parts, and summing up the equations.\n", + "\n", + "It is convenient to think of the PDE system as a vector of equations. The test functions are collected in a vector too, and the variational formulation is the inner product of the vector PDE and the vector test function.\n", + "\n", + "Discretization in time: backward Euler method => time derivative $\\frac{u_i^{n + 1} - u_i^{n}}{\\Delta t}$.\n", + "\n", + "Let $v_1, v_2, v_3$ be teh test functions.\n", + "\n", + "The inner product results in:\n", + "$$\n", + "\\int_\\Omega (\\frac{u_1^{n + 1} - u_1^{n}}{\\Delta t} v_1 + w \\cdot \\nabla u_1^{n+1} v_1 + \\epsilon \\nabla u_1^{n+1} \\cdot \\nabla v_1) dx +\n", + "\\int_\\Omega (\\frac{u_2^{n + 1} - u_2^{n}}{\\Delta t} v_2 + w \\cdot \\nabla u_2^{n+1} v_2 + \\epsilon \\nabla u_2^{n+1} \\cdot \\nabla v_2) dx + \n", + "\\int_\\Omega (\\frac{u_3^{n + 1} - u_3^{n}}{\\Delta t} v_3 + w \\cdot \\nabla u_3^{n+1} v_3 + \\epsilon \\nabla u_3^{n+1} \\cdot \\nabla v_3) dx - \\int_\\Omega (f_1 v_1 + f_2 v_2 + f_3 v_3) dx - \\int_\\Omega (-K u_1^{n+1} u_2^{n+1}v_1 -K u_1^{n+1} u_2^{n+1} v_2 + K u_1^{n+1} u_2^{n+1} v_3 -K u_3^{n+1} v_3) dx = 0\n", + "$$\n", + "\n", + "For this problem it is natural to assume homogeneous Neumann Boundary Conditions on the entire boundary $\\frac{\\partial u_i}{\\partial n} = 0$. The boundary terms vanish when we integrate by parts." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5b14fefa-b413-448a-9a94-14184b51d98a", + "metadata": {}, + "outputs": [], + "source": [ + "import gmsh\n", + "import os\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from tqdm.notebook import tqdm\n", + "from mpi4py import MPI\n", + "from petsc4py import PETSc\n", + "from dolfinx import fem, plot, graph, geometry, io, log\n", + "from dolfinx.fem import petsc as fem_petsc\n", + "from dolfinx.nls import petsc as nls_petsc\n", + "import dolfinx.cpp.mesh as cpp_mesh\n", + "import dolfinx.mesh as dolfin_mesh\n", + "import ufl\n", + "import adios4dolfinx\n", + "import pyvista\n", + "from pyvista.plotting.utilities import active_scalars_algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "dbf5c552-e100-4314-af73-4c7f40b8c511", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Info : Reading 'mesh_flow.msh'...\n", + "Info : 11 entities\n", + "Info : 7584 nodes\n", + "Info : 2064 elements\n", + "Info : Done reading 'mesh_flow.msh'\n" + ] + } + ], + "source": [ + "# Specify the domain (mesh)\n", + "gdim = 2\n", + "gmsh_model_rank = 0\n", + "mesh_comm = MPI.COMM_WORLD\n", + "msh, cell_markers, facet_markers = io.gmshio.read_from_msh(\"mesh_flow.msh\", mesh_comm, gmsh_model_rank, gdim=gdim)\n", + "facet_markers.name = \"Facet markers\"" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "118c4735-e0d2-40f5-b6bc-9f5603450a07", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "be3c4449ee9f47bab2bd88d87b36e328", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Solving Chemistry PDE: 0%| | 0/12800 [00:00