Skip to content

Add advdiffreac.ipynb for demo #181

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add advdiffreac.ipynb for demo #181

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
337 changes: 337 additions & 0 deletions chapter2/advdiffreac.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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<?, ?it/s]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = 0\n",
"T = 8 # Final time\n",
"dt = 1 / 1600 # Time step size\n",
"num_steps = int(T / dt)\n",
"Delta_t = fem.Constant(msh, PETSc.ScalarType(dt))\n",
"eps = fem.Constant(msh, PETSc.ScalarType(0.01)) # diffusivity\n",
"K = fem.Constant(msh, PETSc.ScalarType(10.0)) # rate of reaction\n",
"\n",
"# Choose type of function space\n",
"W = fem.functionspace(\n",
" mesh=msh, \n",
" element=ufl.VectorElement( # For Vector Function Space using Vector Element\n",
" family=\"Lagrange\", \n",
" cell=msh.ufl_cell(), \n",
" degree=2\n",
" )\n",
")\n",
"\n",
"P = ufl.FiniteElement(\n",
" family=\"Lagrange\", \n",
" cell=msh.ufl_cell(), \n",
" degree=1\n",
")\n",
"V = fem.functionspace(\n",
" mesh=msh, \n",
" element=ufl.MixedElement([P, P, P])\n",
")\n",
"\n",
"# Variational Formulation\n",
"# Note that as the problem is non-linear, we have replace the TrialFunction with a Function, \n",
"# which serves as the unknown of our problem\n",
"w = fem.Function(W)\n",
"adios4dolfinx.read_function(w, \"./Flows/flow_function_0\")\n",
"u, (v1, v2, v3) = fem.Function(V), ufl.TestFunctions(V)\n",
"u.x.array[:] = 0.0\n",
"u1, u2, u3 = ufl.split(u)\n",
"\n",
"u_n = fem.Function(V)\n",
"u_n.x.array[:] = 0.0\n",
"u_n1, u_n2, u_n3 = ufl.split(u_n)\n",
"\n",
"f = fem.Function(V)\n",
"f1, f2, f3 = f.sub(0).collapse(), f.sub(1).collapse(), f.sub(2).collapse()\n",
"\n",
"x = ufl.SpatialCoordinate(msh)\n",
"def source(x, target, r):\n",
" result = np.zeros_like(x[0])\n",
" values = (x[0] - target[0])**2 + (x[1] - target[1])**2\n",
" true_mask = values < (r * r)\n",
" false_mask = np.logical_not(true_mask)\n",
" result[true_mask] = 0.1\n",
" result[false_mask] = 0.0\n",
" return result\n",
"f1.interpolate(lambda x : source(x, np.array([0.1, 0.1]), 0.05))\n",
"f2.interpolate(lambda x : source(x, np.array([0.1, 0.3]), 0.05))\n",
"f3.x.array[:] = 0.0\n",
"\n",
"F = ((u1 - u_n1) / Delta_t * v1 + ufl.dot(w, ufl.grad(u1)) * v1 + eps * ufl.dot(ufl.grad(u1), ufl.grad(v1))) * ufl.dx\n",
"F += ((u2 - u_n2) / Delta_t * v2 + ufl.dot(w, ufl.grad(u2)) * v2 + eps * ufl.dot(ufl.grad(u2), ufl.grad(v2))) * ufl.dx\n",
"F += ((u3 - u_n3) / Delta_t * v3 + ufl.dot(w, ufl.grad(u3)) * v3 + eps * ufl.dot(ufl.grad(u3), ufl.grad(v3))) * ufl.dx\n",
"F -= (f1 * v1 + f2 * v2 + f3 * v3) * ufl.dx\n",
"F -= (-K * u1 * u2 * v1 - K * u1 * u2 * v2 + K * u1 * u2 * v3 - K * u3 * v3) * ufl.dx\n",
"\n",
"problem = fem_petsc.NonlinearProblem(F, u)\n",
"\n",
"# Get a solver\n",
"solver = nls_petsc.NewtonSolver(MPI.COMM_WORLD, problem)\n",
"solver.convergence_criterion = \"incremental\"\n",
"solver.rtol = 1e-6\n",
"solver.report = True\n",
"\n",
"ksp = solver.krylov_solver\n",
"opts = PETSc.Options()\n",
"option_prefix = ksp.getOptionsPrefix()\n",
"opts[f\"{option_prefix}ksp_type\"] = \"preonly\"\n",
"opts[f\"{option_prefix}pc_type\"] = \"lu\"\n",
"opts[f\"{option_prefix}pc_factor_mat_solver_type\"] = \"mumps\"\n",
"ksp.setFromOptions()\n",
"\n",
"# log.set_log_level(log.LogLevel.INFO)\n",
"\n",
"# Visualization for gif\n",
"pyvista.start_xvfb()\n",
"V1, V1_to_V = V.sub(0).collapse()\n",
"V2, V2_to_V = V.sub(1).collapse()\n",
"V3, V3_to_V = V.sub(2).collapse()\n",
"\n",
"grid1 = pyvista.UnstructuredGrid(*plot.vtk_mesh(V1))\n",
"grid2 = pyvista.UnstructuredGrid(*plot.vtk_mesh(V2))\n",
"grid3 = pyvista.UnstructuredGrid(*plot.vtk_mesh(V3))\n",
"\n",
"grid1[\"u1\"] = u.sub(0).collapse().x.array\n",
"grid2[\"u2\"] = u.sub(1).collapse().x.array\n",
"grid3[\"u3\"] = u.sub(2).collapse().x.array\n",
"\n",
"plotter = pyvista.Plotter(shape=(3, 1), notebook=True)\n",
"plotter.open_movie(\"u_time_chemistry.mp4\")\n",
"plotter.subplot(0, 0)\n",
"plotter.add_text(f\"time: {t:.2f}\", font_size=12, name=\"timelabel\")\n",
"plotter.add_text(\"Chemical A\", font_size=9, position=\"left_edge\")\n",
"magma = pyvista.LookupTable(cmap=\"magma\")\n",
"magma.scalar_range = (0, 0.03)\n",
"magma.above_range_color = 'r'\n",
"plotter.add_mesh(active_scalars_algorithm(grid1, \"u1\"), show_edges=True, show_scalar_bar=True, cmap=magma)\n",
"plotter.subplot(1, 0)\n",
"plotter.add_text(\"Chemical B\", font_size=9, position=\"left_edge\")\n",
"cool = pyvista.LookupTable(cmap=\"cool\")\n",
"cool.scalar_range = (0, 0.03)\n",
"cool.above_range_color = 'b'\n",
"plotter.add_mesh(active_scalars_algorithm(grid2, \"u2\"), show_edges=True, show_scalar_bar=True, cmap=cool)\n",
"plotter.subplot(2, 0)\n",
"plotter.add_text(\"Chemical C\", font_size=9, position=\"left_edge\")\n",
"lut = pyvista.LookupTable()\n",
"lut.scalar_range = (0, 1.8e-4)\n",
"lut.above_range_color = 'g'\n",
"plotter.add_mesh(active_scalars_algorithm(grid3, \"u3\"), show_edges=True, show_scalar_bar=True, cmap=lut)\n",
"plotter.link_views()\n",
"plotter.view_xy()\n",
"plotter.camera.zoom(3)\n",
"\n",
"progress = tqdm(desc=\"Solving Chemistry PDE\", total=num_steps)\n",
"for i in range(num_steps):\n",
" progress.update(1)\n",
" adios4dolfinx.read_function(w, f\"./Flows/flow_function_{i}\")\n",
"\n",
" try:\n",
" n, converged = solver.solve(u)\n",
" except Exception as error:\n",
" print(\"An error occurred:\", type(error).__name__) # An error occurred: NameError\n",
" continue\n",
" # print(f\"Step {i+1}: num iterations: {n}\")\n",
"\n",
" # Update variable with solution form this time step\n",
" u_n.x.array[:] = u.x.array\n",
" u.x.scatter_forward()\n",
" # For animation\n",
" grid1[\"u1\"] = u.sub(0).collapse().x.array\n",
" grid2[\"u2\"] = u.sub(1).collapse().x.array\n",
" grid3[\"u3\"] = u.sub(2).collapse().x.array\n",
" plotter.subplot(0, 0)\n",
" plotter.add_text(f\"time: {t:.2f}\", font_size=12, name=\"timelabel\")\n",
" plotter.render()\n",
" plotter.subplot(1, 0)\n",
" plotter.render()\n",
" plotter.subplot(2, 0)\n",
" plotter.render()\n",
" plotter.write_frame()\n",
" t += dt\n",
" \n",
"plotter.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a484f495-3ae7-4980-a5ec-31a769a2c622",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.8.18"
}
},
"nbformat": 4,
"nbformat_minor": 5
}