7-TransientNR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "21cf672e-ba27-480c-a1e0-39c6dddd09f8",
   "metadata": {},
   "source": [
    "# Your first nonlinear transient solver- nonlinear diffusion\n",
    "\n",
    "- [x] newton-raphson method\n",
    "- [x] weak form\n",
    "- [x] residual and jacobian\n",
    "- [x] nonlinear iteration\n",
    "- [x] time stepping\n",
    "- [x] convergence criterion\n",
    "- [x] nonlinear diffusion equation\n",
    "\n",
    "Author: Yang Bai @ M3 Group\n",
    "\n",
    "Date  : 2022.06.04\n",
    "\n",
    "QQ group: 628204857"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "279b0813-25ab-4118-b26b-f08ef9a516e0",
   "metadata": {},
   "source": [
    "# Governing equation\n",
    "\n",
    "$$\\frac{\\partial c}{\\partial t}=\\nabla\\cdot(D(c)\\nabla c)+f(c)$$, with $-D\\nabla c\\cdot\\vec{n}=j_{0}$\n",
    "\n",
    "## Weak form\n",
    "\n",
    "$$\\int_{\\Omega}\\dot{c}\\delta c dV=\\int_{\\partial\\Omega}D(c)\\nabla c\\cdot\\vec{n}\\delta c dS-\\int_{\\Omega}D\\nabla c\\nabla\\delta cdV+\\int_{\\Omega}f(c)\\delta c dV$$\n",
    "\n",
    "### Residual\n",
    "$$R_{c}^{I}=\\int_{\\Omega}\\dot{c}N^{I}dV+\\int_{\\Omega}D(c)\\nabla c\\nabla N^{I}dV-\\int_{\\Omega}f(c)N^{I}dV+\\int_{\\partial\\Omega}j_{0}N^{I}dS$$\n",
    "\n",
    "### Jacobian\n",
    "$$\\begin{aligned}\n",
    "K_{cc}^{IJ}&=-\\frac{\\partial R_{c}^{I}}{\\partial c^{J}}\\\\\n",
    "&=-\\int_{\\Omega}\\frac{\\partial\\dot{c}}{\\partial c}N^{J}N^{I}dV\\\\\n",
    "&-\\int_{\\Omega}\\frac{\\partial D(c)}{\\partial c}N^{J}\\nabla c\\nabla N^{I} dV\n",
    "-\\int_{\\Omega}D(c)\\nabla N^{J}\\nabla N^{I}dV\\\\\n",
    "&+\\int_{\\Omega}\\frac{\\partial f(c)}{\\partial c}N^{J}N^{I}dV\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f94144e7-81f3-4fc6-b65f-21395e167413",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "write result to nldiff-000000.vtu\n",
      "step=     1, dt=   2.00000e-03, total time=   2.00000e-03\n",
      "    iters=  1, |R0 |=   1.13069e-01, |R |=   1.13069e-01\n",
      "    iters=  2, |R0 |=   1.13069e-01, |R |=   2.12859e-05\n",
      "    iters=  3, |R0 |=   1.13069e-01, |R |=   2.16693e-11\n",
      "step=     2, dt=   2.00000e-03, total time=   4.00000e-03\n",
      "    iters=  1, |R0 |=   1.12965e-01, |R |=   1.12965e-01\n",
      "    iters=  2, |R0 |=   1.12965e-01, |R |=   7.48628e-06\n",
      "    iters=  3, |R0 |=   1.12965e-01, |R |=   3.05335e-12\n",
      "step=     3, dt=   2.00000e-03, total time=   6.00000e-03\n",
      "    iters=  1, |R0 |=   1.12959e-01, |R |=   1.12959e-01\n",
      "    iters=  2, |R0 |=   1.12959e-01, |R |=   5.44658e-06\n",
      "    iters=  3, |R0 |=   1.12959e-01, |R |=   5.76423e-13\n",
      "step=     4, dt=   2.00000e-03, total time=   8.00000e-03\n",
      "    iters=  1, |R0 |=   1.12971e-01, |R |=   1.12971e-01\n",
      "    iters=  2, |R0 |=   1.12971e-01, |R |=   4.95716e-06\n",
      "    iters=  3, |R0 |=   1.12971e-01, |R |=   1.79581e-13\n",
      "step=     5, dt=   2.00000e-03, total time=   1.00000e-02\n",
      "    iters=  1, |R0 |=   1.12995e-01, |R |=   1.12995e-01\n",
      "    iters=  2, |R0 |=   1.12995e-01, |R |=   4.76764e-06\n",
      "    iters=  3, |R0 |=   1.12995e-01, |R |=   9.24800e-14\n",
      "write result to nldiff-000005.vtu\n",
      "---------------------------------------------------------\n",
      "step=     6, dt=   2.00000e-03, total time=   1.20000e-02\n",
      "    iters=  1, |R0 |=   1.13031e-01, |R |=   1.13031e-01\n",
      "    iters=  2, |R0 |=   1.13031e-01, |R |=   4.67827e-06\n",
      "    iters=  3, |R0 |=   1.13031e-01, |R |=   6.02366e-14\n",
      "step=     7, dt=   2.00000e-03, total time=   1.40000e-02\n",
      "    iters=  1, |R0 |=   1.13077e-01, |R |=   1.13077e-01\n",
      "    iters=  2, |R0 |=   1.13077e-01, |R |=   4.63198e-06\n",
      "    iters=  3, |R0 |=   1.13077e-01, |R |=   4.43947e-14\n",
      "step=     8, dt=   2.00000e-03, total time=   1.60000e-02\n",
      "    iters=  1, |R0 |=   1.13133e-01, |R |=   1.13133e-01\n",
      "    iters=  2, |R0 |=   1.13133e-01, |R |=   4.60710e-06\n",
      "    iters=  3, |R0 |=   1.13133e-01, |R |=   3.56324e-14\n",
      "step=     9, dt=   2.00000e-03, total time=   1.80000e-02\n",
      "    iters=  1, |R0 |=   1.13199e-01, |R |=   1.13199e-01\n",
      "    iters=  2, |R0 |=   1.13199e-01, |R |=   4.59410e-06\n",
      "    iters=  3, |R0 |=   1.13199e-01, |R |=   3.04532e-14\n",
      "step=    10, dt=   2.00000e-03, total time=   2.00000e-02\n",
      "    iters=  1, |R0 |=   1.13274e-01, |R |=   1.13274e-01\n",
      "    iters=  2, |R0 |=   1.13274e-01, |R |=   4.58830e-06\n",
      "    iters=  3, |R0 |=   1.13274e-01, |R |=   2.72252e-14\n",
      "write result to nldiff-000010.vtu\n",
      "---------------------------------------------------------\n",
      "step=    11, dt=   2.00000e-03, total time=   2.20000e-02\n",
      "    iters=  1, |R0 |=   1.13359e-01, |R |=   1.13359e-01\n",
      "    iters=  2, |R0 |=   1.13359e-01, |R |=   4.58721e-06\n",
      "    iters=  3, |R0 |=   1.13359e-01, |R |=   2.51231e-14\n",
      "step=    12, dt=   2.00000e-03, total time=   2.40000e-02\n",
      "    iters=  1, |R0 |=   1.13453e-01, |R |=   1.13453e-01\n",
      "    iters=  2, |R0 |=   1.13453e-01, |R |=   4.58943e-06\n",
      "    iters=  3, |R0 |=   1.13453e-01, |R |=   2.37026e-14\n",
      "step=    13, dt=   2.00000e-03, total time=   2.60000e-02\n",
      "    iters=  1, |R0 |=   1.13557e-01, |R |=   1.13557e-01\n",
      "    iters=  2, |R0 |=   1.13557e-01, |R |=   4.59410e-06\n",
      "    iters=  3, |R0 |=   1.13557e-01, |R |=   2.27100e-14\n",
      "step=    14, dt=   2.00000e-03, total time=   2.80000e-02\n",
      "    iters=  1, |R0 |=   1.13670e-01, |R |=   1.13670e-01\n",
      "    iters=  2, |R0 |=   1.13670e-01, |R |=   4.60071e-06\n",
      "    iters=  3, |R0 |=   1.13670e-01, |R |=   2.19831e-14\n",
      "step=    15, dt=   2.00000e-03, total time=   3.00000e-02\n",
      "    iters=  1, |R0 |=   1.13792e-01, |R |=   1.13792e-01\n",
      "    iters=  2, |R0 |=   1.13792e-01, |R |=   4.60891e-06\n",
      "    iters=  3, |R0 |=   1.13792e-01, |R |=   2.14399e-14\n",
      "write result to nldiff-000015.vtu\n",
      "---------------------------------------------------------\n",
      "step=    16, dt=   2.00000e-03, total time=   3.20000e-02\n",
      "    iters=  1, |R0 |=   1.13924e-01, |R |=   1.13924e-01\n",
      "    iters=  2, |R0 |=   1.13924e-01, |R |=   4.61847e-06\n",
      "    iters=  3, |R0 |=   1.13924e-01, |R |=   2.10235e-14\n",
      "step=    17, dt=   2.00000e-03, total time=   3.40000e-02\n",
      "    iters=  1, |R0 |=   1.14066e-01, |R |=   1.14066e-01\n",
      "    iters=  2, |R0 |=   1.14066e-01, |R |=   4.62925e-06\n",
      "    iters=  3, |R0 |=   1.14066e-01, |R |=   2.06974e-14\n",
      "step=    18, dt=   2.00000e-03, total time=   3.60000e-02\n",
      "    iters=  1, |R0 |=   1.14217e-01, |R |=   1.14217e-01\n",
      "    iters=  2, |R0 |=   1.14217e-01, |R |=   4.64113e-06\n",
      "    iters=  3, |R0 |=   1.14217e-01, |R |=   2.04409e-14\n",
      "step=    19, dt=   2.00000e-03, total time=   3.80000e-02\n",
      "    iters=  1, |R0 |=   1.14377e-01, |R |=   1.14377e-01\n",
      "    iters=  2, |R0 |=   1.14377e-01, |R |=   4.65406e-06\n",
      "    iters=  3, |R0 |=   1.14377e-01, |R |=   2.02405e-14\n",
      "step=    20, dt=   2.00000e-03, total time=   4.00000e-02\n",
      "    iters=  1, |R0 |=   1.14547e-01, |R |=   1.14547e-01\n",
      "    iters=  2, |R0 |=   1.14547e-01, |R |=   4.66797e-06\n",
      "    iters=  3, |R0 |=   1.14547e-01, |R |=   2.00634e-14\n",
      "write result to nldiff-000020.vtu\n",
      "---------------------------------------------------------\n",
      "step=    21, dt=   2.00000e-03, total time=   4.20000e-02\n",
      "    iters=  1, |R0 |=   1.14727e-01, |R |=   1.14727e-01\n",
      "    iters=  2, |R0 |=   1.14727e-01, |R |=   4.68283e-06\n",
      "    iters=  3, |R0 |=   1.14727e-01, |R |=   1.99043e-14\n",
      "step=    22, dt=   2.00000e-03, total time=   4.40000e-02\n",
      "    iters=  1, |R0 |=   1.14916e-01, |R |=   1.14916e-01\n",
      "    iters=  2, |R0 |=   1.14916e-01, |R |=   4.69862e-06\n",
      "    iters=  3, |R0 |=   1.14916e-01, |R |=   1.97837e-14\n",
      "step=    23, dt=   2.00000e-03, total time=   4.60000e-02\n",
      "    iters=  1, |R0 |=   1.15116e-01, |R |=   1.15116e-01\n",
      "    iters=  2, |R0 |=   1.15116e-01, |R |=   4.71532e-06\n",
      "    iters=  3, |R0 |=   1.15116e-01, |R |=   1.96862e-14\n",
      "step=    24, dt=   2.00000e-03, total time=   4.80000e-02\n",
      "    iters=  1, |R0 |=   1.15325e-01, |R |=   1.15325e-01\n",
      "    iters=  2, |R0 |=   1.15325e-01, |R |=   4.73294e-06\n",
      "    iters=  3, |R0 |=   1.15325e-01, |R |=   1.96174e-14\n",
      "step=    25, dt=   2.00000e-03, total time=   5.00000e-02\n",
      "    iters=  1, |R0 |=   1.15543e-01, |R |=   1.15543e-01\n",
      "    iters=  2, |R0 |=   1.15543e-01, |R |=   4.75145e-06\n",
      "    iters=  3, |R0 |=   1.15543e-01, |R |=   1.95382e-14\n",
      "write result to nldiff-000025.vtu\n",
      "---------------------------------------------------------\n",
      "step=    26, dt=   2.00000e-03, total time=   5.20000e-02\n",
      "    iters=  1, |R0 |=   1.15772e-01, |R |=   1.15772e-01\n",
      "    iters=  2, |R0 |=   1.15772e-01, |R |=   4.77086e-06\n",
      "    iters=  3, |R0 |=   1.15772e-01, |R |=   1.95029e-14\n",
      "step=    27, dt=   2.00000e-03, total time=   5.40000e-02\n",
      "    iters=  1, |R0 |=   1.16011e-01, |R |=   1.16011e-01\n",
      "    iters=  2, |R0 |=   1.16011e-01, |R |=   4.79118e-06\n",
      "    iters=  3, |R0 |=   1.16011e-01, |R |=   1.94855e-14\n",
      "step=    28, dt=   2.00000e-03, total time=   5.60000e-02\n",
      "    iters=  1, |R0 |=   1.16260e-01, |R |=   1.16260e-01\n",
      "    iters=  2, |R0 |=   1.16260e-01, |R |=   4.81242e-06\n",
      "    iters=  3, |R0 |=   1.16260e-01, |R |=   1.94438e-14\n",
      "step=    29, dt=   2.00000e-03, total time=   5.80000e-02\n",
      "    iters=  1, |R0 |=   1.16519e-01, |R |=   1.16519e-01\n",
      "    iters=  2, |R0 |=   1.16519e-01, |R |=   4.83457e-06\n",
      "    iters=  3, |R0 |=   1.16519e-01, |R |=   1.93983e-14\n",
      "step=    30, dt=   2.00000e-03, total time=   6.00000e-02\n",
      "    iters=  1, |R0 |=   1.16788e-01, |R |=   1.16788e-01\n",
      "    iters=  2, |R0 |=   1.16788e-01, |R |=   4.85764e-06\n",
      "    iters=  3, |R0 |=   1.16788e-01, |R |=   1.93792e-14\n",
      "write result to nldiff-000030.vtu\n",
      "---------------------------------------------------------\n",
      "step=    31, dt=   2.00000e-03, total time=   6.20000e-02\n",
      "    iters=  1, |R0 |=   1.17068e-01, |R |=   1.17068e-01\n",
      "    iters=  2, |R0 |=   1.17068e-01, |R |=   4.88166e-06\n",
      "    iters=  3, |R0 |=   1.17068e-01, |R |=   1.93926e-14\n",
      "step=    32, dt=   2.00000e-03, total time=   6.40000e-02\n",
      "    iters=  1, |R0 |=   1.17358e-01, |R |=   1.17358e-01\n",
      "    iters=  2, |R0 |=   1.17358e-01, |R |=   4.90663e-06\n",
      "    iters=  3, |R0 |=   1.17358e-01, |R |=   1.94064e-14\n",
      "step=    33, dt=   2.00000e-03, total time=   6.60000e-02\n",
      "    iters=  1, |R0 |=   1.17658e-01, |R |=   1.17658e-01\n",
      "    iters=  2, |R0 |=   1.17658e-01, |R |=   4.93256e-06\n",
      "    iters=  3, |R0 |=   1.17658e-01, |R |=   1.94326e-14\n",
      "step=    34, dt=   2.00000e-03, total time=   6.80000e-02\n",
      "    iters=  1, |R0 |=   1.17970e-01, |R |=   1.17970e-01\n",
      "    iters=  2, |R0 |=   1.17970e-01, |R |=   4.95947e-06\n",
      "    iters=  3, |R0 |=   1.17970e-01, |R |=   1.94339e-14\n",
      "step=    35, dt=   2.00000e-03, total time=   7.00000e-02\n",
      "    iters=  1, |R0 |=   1.18292e-01, |R |=   1.18292e-01\n",
      "    iters=  2, |R0 |=   1.18292e-01, |R |=   4.98737e-06\n",
      "    iters=  3, |R0 |=   1.18292e-01, |R |=   1.94976e-14\n",
      "write result to nldiff-000035.vtu\n",
      "---------------------------------------------------------\n",
      "step=    36, dt=   2.00000e-03, total time=   7.20000e-02\n",
      "    iters=  1, |R0 |=   1.18625e-01, |R |=   1.18625e-01\n",
      "    iters=  2, |R0 |=   1.18625e-01, |R |=   5.01628e-06\n",
      "    iters=  3, |R0 |=   1.18625e-01, |R |=   1.95438e-14\n",
      "step=    37, dt=   2.00000e-03, total time=   7.40000e-02\n",
      "    iters=  1, |R0 |=   1.18969e-01, |R |=   1.18969e-01\n",
      "    iters=  2, |R0 |=   1.18969e-01, |R |=   5.04622e-06\n",
      "    iters=  3, |R0 |=   1.18969e-01, |R |=   1.95886e-14\n",
      "step=    38, dt=   2.00000e-03, total time=   7.60000e-02\n",
      "    iters=  1, |R0 |=   1.19324e-01, |R |=   1.19324e-01\n",
      "    iters=  2, |R0 |=   1.19324e-01, |R |=   5.07721e-06\n",
      "    iters=  3, |R0 |=   1.19324e-01, |R |=   1.96322e-14\n",
      "step=    39, dt=   2.00000e-03, total time=   7.80000e-02\n",
      "    iters=  1, |R0 |=   1.19690e-01, |R |=   1.19690e-01\n",
      "    iters=  2, |R0 |=   1.19690e-01, |R |=   5.10926e-06\n",
      "    iters=  3, |R0 |=   1.19690e-01, |R |=   1.97176e-14\n",
      "step=    40, dt=   2.00000e-03, total time=   8.00000e-02\n",
      "    iters=  1, |R0 |=   1.20067e-01, |R |=   1.20067e-01\n",
      "    iters=  2, |R0 |=   1.20067e-01, |R |=   5.14239e-06\n",
      "    iters=  3, |R0 |=   1.20067e-01, |R |=   1.97652e-14\n",
      "write result to nldiff-000040.vtu\n",
      "---------------------------------------------------------\n",
      "step=    41, dt=   2.00000e-03, total time=   8.20000e-02\n",
      "    iters=  1, |R0 |=   1.20456e-01, |R |=   1.20456e-01\n",
      "    iters=  2, |R0 |=   1.20456e-01, |R |=   5.17663e-06\n",
      "    iters=  3, |R0 |=   1.20456e-01, |R |=   1.98688e-14\n",
      "step=    42, dt=   2.00000e-03, total time=   8.40000e-02\n",
      "    iters=  1, |R0 |=   1.20857e-01, |R |=   1.20857e-01\n",
      "    iters=  2, |R0 |=   1.20857e-01, |R |=   5.21200e-06\n",
      "    iters=  3, |R0 |=   1.20857e-01, |R |=   1.99779e-14\n",
      "step=    43, dt=   2.00000e-03, total time=   8.60000e-02\n",
      "    iters=  1, |R0 |=   1.21270e-01, |R |=   1.21270e-01\n",
      "    iters=  2, |R0 |=   1.21270e-01, |R |=   5.24851e-06\n",
      "    iters=  3, |R0 |=   1.21270e-01, |R |=   2.00173e-14\n",
      "step=    44, dt=   2.00000e-03, total time=   8.80000e-02\n",
      "    iters=  1, |R0 |=   1.21694e-01, |R |=   1.21694e-01\n",
      "    iters=  2, |R0 |=   1.21694e-01, |R |=   5.28620e-06\n",
      "    iters=  3, |R0 |=   1.21694e-01, |R |=   2.01637e-14\n",
      "step=    45, dt=   2.00000e-03, total time=   9.00000e-02\n",
      "    iters=  1, |R0 |=   1.22130e-01, |R |=   1.22130e-01\n",
      "    iters=  2, |R0 |=   1.22130e-01, |R |=   5.32509e-06\n",
      "    iters=  3, |R0 |=   1.22130e-01, |R |=   2.02108e-14\n",
      "write result to nldiff-000045.vtu\n",
      "---------------------------------------------------------\n",
      "step=    46, dt=   2.00000e-03, total time=   9.20000e-02\n",
      "    iters=  1, |R0 |=   1.22579e-01, |R |=   1.22579e-01\n",
      "    iters=  2, |R0 |=   1.22579e-01, |R |=   5.36520e-06\n",
      "    iters=  3, |R0 |=   1.22579e-01, |R |=   2.03732e-14\n",
      "step=    47, dt=   2.00000e-03, total time=   9.40000e-02\n",
      "    iters=  1, |R0 |=   1.23040e-01, |R |=   1.23040e-01\n",
      "    iters=  2, |R0 |=   1.23040e-01, |R |=   5.40657e-06\n",
      "    iters=  3, |R0 |=   1.23040e-01, |R |=   2.04949e-14\n",
      "step=    48, dt=   2.00000e-03, total time=   9.60000e-02\n",
      "    iters=  1, |R0 |=   1.23514e-01, |R |=   1.23514e-01\n",
      "    iters=  2, |R0 |=   1.23514e-01, |R |=   5.44921e-06\n",
      "    iters=  3, |R0 |=   1.23514e-01, |R |=   2.06845e-14\n",
      "step=    49, dt=   2.00000e-03, total time=   9.80000e-02\n",
      "    iters=  1, |R0 |=   1.24000e-01, |R |=   1.24000e-01\n",
      "    iters=  2, |R0 |=   1.24000e-01, |R |=   5.49316e-06\n",
      "    iters=  3, |R0 |=   1.24000e-01, |R |=   2.07789e-14\n",
      "write result to nldiff-000049.vtu\n",
      "---------------------------------------------------------\n",
      "simulation is done!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARqElEQVR4nO3dbYxcZ3nG8f+FSQQIVLc1pZbt1G7qvrgQ3LCk6RtxgbZ2PtStxAfHVSICyAoiFVVVQdQPoIpKFdA3oQYsi1oQqcFCIqVu6jRCtEmKaMAONSFOmnQxgmwdKQ4oCYG+yHD3w4ybedZre+PMnj0z+/9JK+2cc3bv++xDfHHmzNyTqkKSpNNesNwNSJL6xWCQJDUMBklSw2CQJDUMBklS44XL3cDztWbNmtq4ceNytyFJE+W+++57oqpevtC+iQ+GjRs3cuTIkeVuQ5ImSpKvn22fTyVJkhoGgySpYTBIkhoGgySpYTBIkhqdBUOS/UkeT/LAWfYnyYeSzCa5P8nlXfUmSXpWl1cMHwO2n2P/DmDz8GsP8JEOepIkzdPZ+xiq6p4kG89xyE7glhrMAb83yeoka6vqsaXoZ9++fdx6660L7jtx4gRJWLt27Vj2jfv3TUIf1pqePqzVz1oAu3fvZs+ePQvuez769Aa3dcCjI4/nhtvOCIYkexhcVXDJJZdcULFbb72Vo0ePsnXr1jP2HT9+HGDBxbiQfeP+fZPQh7Wmpw9r9bPW0aNHAZYkGNLlB/UMrxhur6pXLrDvH4A/qarPDR9/FnhXVd13rt85MzNTF/LO523btgFw1113nbFv9erVADz55JNj2Tfu3zcJfVhrevqwVj9rnevfsMVIcl9VzSy0r0+vSpoDNow8Xg+cWKZeJGnF6lMwHASuG7466UrgqaW6vyBJOrvO7jEk+QSwDViTZA54L3ARQFXtBQ4BVwOzwHeB67vqTZL0rC5flXTNefYX8I6O2pEknUWfnkqSJPWAwSBJahgMkqSGwSBJahgMkqSGwSBJavRpVlKn7r77buDZt5yPeuqpp8a6b9y/bxL6sNb09GGtftZ65plnuPTSS8/YPg5eMUjShFqqWXcr9orhqquuAhyiZ63+1upLH9bqZ63TQ/SWglcMkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqSGwSBJahgMkqTGip2V5HRVa/W9Vl/6sFY/azldVZJ0BqerjpnTVa3V91p96cNa/azldFVJUmcMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSo9NgSLI9ycNJZpPctMD+H0jy90m+nORYkuu77E+S1GEwJFkF3AzsALYA1yTZMu+wdwAPVtWrgW3AnyW5uKseJUndjsS4ApitquMASQ4AO4EHR44p4GVJArwU+BZwaimacYietfpeqy99WKuftaZliN464NGRx3PDbaP+CvgZ4ATwFeCdVfX9+b8oyZ4kR5IcOXny5FL1K0m9Ng1D9LLAtvln9RvAUeD1wKXAZ5L8S1U93fxQ1T5gH8DMzMwF/WUcometvtfqSx/W6metaRmiNwdsGHm8nsGVwajrgdtqYBb4GvDTHfUnSaLbYDgMbE6yaXhDeRdwcN4x3wDeAJDkFcBPAcc77FGSVrzOnkqqqlNJbgTuBFYB+6vqWJIbhvv3Au8DPpbkKwyeenp3VT3RVY+SpI4/qKeqDgGH5m3bO/L9CeDXu+xJktTync+SpIbBIElqGAySpIbBIElqGAySpIbBIElqdPpy1T5xiJ61+l6rL31Yq5+1pmWIniRpjKZhiF6vOETPWn2v1Zc+rNXPWtMyRE+SNAEMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDUMBklSY8WOxHCInrX6XqsvfVirn7UcoidJOoND9MbMIXrW6nutvvRhrX7WcoieJKkzBoMkqWEwSJIaBoMkqWEwSJIaBoMkqWEwSJIaBoMkqdFpMCTZnuThJLNJbjrLMduSHE1yLMndXfYnSerwnc9JVgE3A78GzAGHkxysqgdHjlkNfBjYXlXfSPIjXfUnSRrociTGFcBsVR0HSHIA2Ak8OHLMbuC2qvoGQFU9vlTNOETPWn2v1Zc+rNXPWtMyRG8d8OjI47nhtlE/CfxgkruS3JfkuoV+UZI9SY4kOXLy5MklaleS+m0ahuhlgW3zz+qFwGuANwAvBv41yb1V9UjzQ1X7gH0AMzMzF/SXcYietfpeqy99WKuftZZyiF6XwTAHbBh5vB44scAxT1TVd4DvJLkHeDXwCJKkTnT5VNJhYHOSTUkuBnYBB+cd83fAryR5YZKXAD8PPNRhj5K04nV2xVBVp5LcCNwJrAL2V9WxJDcM9++tqoeS/CNwP/B94KNV9UBXPUqSOv6gnqo6BByat23vvMcfBD7YZV+SpGf5zmdJUsNgkCQ1DAZJUsNgkCQ1DAZJUsNgkCQ1On25ap84RM9afa/Vlz6s1c9a0zJET5I0RtMwRK9XHKJnrb7X6ksf1upnraUcoucVgySpYTBIkhrnDIYkV3bViCSpH853xXBPkvclWbH3IiRppTlfMOwArgW+mGRLB/1IkpbZOYOhqj4LvAr4N+BIkt/vpCtJ0rI5783nqvp2Vb2VwZXDB5I8k+Tp0a+lb1OS1JVF3TtIMgP8MfAfwJ8Cp5ayKUnS8jlnMAxvOr8XeDfwYeCmqvrvLhqTJC2P810xHAZ+CNgxvN8gSZpy57vH8ABwmaEgSSvHOa8YqurarhrpmtNVrdX3Wn3pw1r9rOV0VUnSGZyuOmZOV7VW32v1pQ9r9bOW01UlSZ0xGCRJDYNBktQwGCRJDYNBktQwGCRJDYNBktQwGCRJjU6DIcn2JA8nmU1y0zmOe22S7yV5U5f9SZI6DIYkq4CbGXxc6BbgmoU+LnR43PuBO7vqTZL0rC5HYlwBzFbVcYAkB4CdwIPzjvtd4FPAa5eyGYfoWavvtfrSh7X6WWtahuitAx4deTw33Pb/kqwDfhvYe65flGRPkiNJjpw8eXLsjUrSJJiGIXpZYNv8s/pL4N1V9b1kocOHP1S1D9gHMDMzc0F/GYfoWavvtfrSh7X6WWsph+h1GQxzwIaRx+uBE/OOmQEODENhDXB1klNV9elOOpQkdRoMh4HNSTYB/wnsAnaPHlBVm05/n+RjwO2GgiR1q7NgqKpTSW5k8GqjVcD+qjqW5Ibh/nPeV5AkdaPTD+qpqkPAoXnbFgyEqnpzFz1Jklq+81mS1DAYJEkNg0GS1DAYJEkNg0GS1DAYJEmNTl+u2icO0bNW32v1pQ9r9bPWtAzRkySN0TQM0esVh+hZq++1+tKHtfpZaymH6HnFIElqGAySpIbBIElqGAySpIbBIElqGAySpIbBIElqGAySpIbBIElqGAySpMaKHYnhED1r9b1WX/qwVj9rOURPknQGh+iNmUP0rNX3Wn3pw1r9rOUQPUlSZwwGSVLDYJAkNQwGSVLDYJAkNQwGSVLDYJAkNQwGSVKj02BIsj3Jw0lmk9y0wP7fSXL/8OvzSV7dZX+SpA6DIckq4GZgB7AFuCbJlnmHfQ24qqouA94H7OuqP0nSQJcjMa4AZqvqOECSA8BO4MHTB1TV50eOvxdYv1TNOETPWn2v1Zc+rNXPWtMyRG8d8OjI47nhtrN5K3DHQjuS7ElyJMmRkydPjrFFSZoc0zBELwtsW/Cskvwqg2D45YX2V9U+hk8zzczMXNBfxiF61up7rb70Ya1+1lrKIXpdBsMcsGHk8XrgxPyDklwGfBTYUVXf7Kg3SdJQl08lHQY2J9mU5GJgF3Bw9IAklwC3AddW1SMd9iZJGursiqGqTiW5EbgTWAXsr6pjSW4Y7t8LvAf4YeDDSQBOVdVMVz1Kkjr+oJ6qOgQcmrdt78j3bwPe1mVPkqSW73yWJDUMBklSw2CQJDUMBklSw2CQJDUMBklSw2CQJDU6fR9Dnzhd1Vp9r9WXPqzVz1rTMl1VkjRG0zBdtVecrmqtvtfqSx/W6metpZyu6hWDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGgaDJKlhMEiSGit2JIZD9KzV91p96cNa/azlED1J0hkcojdmDtGzVt9r9aUPa/WzlkP0JEmdMRgkSQ2DQZLUMBgkSQ2DQZLUMBgkSQ2DQZLUMBgkSY1OgyHJ9iQPJ5lNctMC+5PkQ8P99ye5vMv+JEkdBkOSVcDNwA5gC3BNki3zDtsBbB5+7QE+0lV/kqSBLkdiXAHMVtVxgCQHgJ3AgyPH7ARuqcEAkHuTrE6ytqoeG3czp4foLfS28tODq8a1b9y/bxL6sNb09GGtftY6evQoW7duPWP7OHT5VNI64NGRx3PDbc/1GJLsSXIkyZGTJ0+OvVFJ6rutW7eye/fuJfndXV4xZIFt80cDLuYYqmofsA9gZmbmgsYLLtVUQkmadF1eMcwBG0YerwdOXMAxkqQl1GUwHAY2J9mU5GJgF3Bw3jEHgeuGr066EnhqKe4vSJLOrrOnkqrqVJIbgTuBVcD+qjqW5Ibh/r3AIeBqYBb4LnB9V/1JkgY6/aCeqjrE4B//0W17R74v4B1d9iRJavnOZ0lSw2CQJDUMBklSw2CQJDUy6W/0SnIS+PoF/vga4IkxtrOcPJd+mpZzmZbzAM/ltB+rqpcvtGPig+H5SHKkqmaWu49x8Fz6aVrOZVrOAzyXxfCpJElSw2CQJDVWejDsW+4Gxshz6adpOZdpOQ/wXM5rRd9jkCSdaaVfMUiS5jEYJEmNFREMSbYneTjJbJKbFtifJB8a7r8/yeXL0ediLOJctiV5KsnR4dd7lqPP80myP8njSR44y/5JWpPzncukrMmGJP+c5KEkx5K8c4FjJmJdFnkuk7IuL0ryxSRfHp7LHy1wzHjXpaqm+ovBiO+vAj8OXAx8Gdgy75irgTsYfILclcAXlrvv53Eu24Dbl7vXRZzL64DLgQfOsn8i1mSR5zIpa7IWuHz4/cuARyb4v5XFnMukrEuAlw6/vwj4AnDlUq7LSrhiuAKYrarjVfW/wAFg57xjdgK31MC9wOoka7tudBEWcy4ToaruAb51jkMmZU0Wcy4Toaoeq6ovDb//NvAQZ37m+kSsyyLPZSIM/9bPDB9eNPya/6qhsa7LSgiGdcCjI4/nOPN/IIs5pg8W2+cvDC8770jys920NnaTsiaLNVFrkmQj8HMM/t/pqIlbl3OcC0zIuiRZleQo8Djwmapa0nXp9IN6lkkW2DY/bRdzTB8sps8vMZiB8kySq4FPA5uXurElMClrshgTtSZJXgp8Cvi9qnp6/u4FfqS363Kec5mYdamq7wFbk6wG/jbJK6tq9J7WWNdlJVwxzAEbRh6vB05cwDF9cN4+q+rp05edNfjEvIuSrOmuxbGZlDU5r0lakyQXMfiH9G+q6rYFDpmYdTnfuUzSupxWVU8CdwHb5+0a67qshGA4DGxOsinJxcAu4OC8Yw4C1w3v7F8JPFVVj3Xd6CKc91yS/GiSDL+/gsEaf7PzTp+/SVmT85qUNRn2+NfAQ1X152c5bCLWZTHnMkHr8vLhlQJJXgy8Efj3eYeNdV2m/qmkqjqV5EbgTgav6tlfVceS3DDcv5fB51BfDcwC3wWuX65+z2WR5/Im4O1JTgH/Beyq4csW+iTJJxi8KmRNkjngvQxuqk3UmsCizmUi1gT4JeBa4CvD57MB/hC4BCZuXRZzLpOyLmuBjydZxSC8PllVty/lv2GOxJAkNVbCU0mSpOfAYJAkNQwGSVLDYJAkNQwGSVLDYJAkNQwGaUySvCDJPUnmv+nwJRmMSv/IcvUmPRcGgzQmVfV94M3A65O8ZWTX+xm8mfQPlqMv6bnyDW7SmA3fkfoB4FXATzB4p/q2qvrcsjYmLZLBIC2BJHcCLwY2Ageq6l3L25G0eAaDtASSbGLwaXtfBV5ZVf+zzC1Ji+Y9BmlpvIXBYLb1DD6KVZoYXjFIY5bktcDngd8E3g68AvjF4YetSL3nFYM0RkleBNwCfKyq7gD2MLgB7T0GTQyvGKQxSvIXwG8Blw0/hJ4ku4CPA6+Z93GMUi8ZDNKYJHkd8E/AG6vqrnn7PsngXsOVVXVqGdqTFs1gkCQ1vMcgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkhsEgSWoYDJKkxv8B97LHCYXiHkMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from FEToy.mesh.lagrange2dmesh import mesh2d\n",
    "from FEToy.fe.shapefun import shape2d\n",
    "from FEToy.fe.gaussrule import gausspoint2d\n",
    "### for surface integration\n",
    "from FEToy.mesh.lagrange1dmesh import mesh1d\n",
    "from FEToy.fe.shapefun import shape1d\n",
    "from FEToy.fe.gaussrule import gausspoint1d\n",
    "### for result output\n",
    "from FEToy.postprocess.Result import ResultIO\n",
    "import numpy as np\n",
    "\n",
    "# define your mesh\n",
    "mymesh=mesh2d(xmax=3.0,ymax=1.0,nx=45,ny=15,meshtype='quad4')\n",
    "mymesh.createmesh()\n",
    "mymesh.plotmesh()\n",
    "\n",
    "# define your shape function\n",
    "shp=shape2d(meshtype='quad4')\n",
    "shp.update()\n",
    "# for surface integration\n",
    "shp1d=shape1d(meshtype='edge2')\n",
    "shp1d.update()\n",
    "\n",
    "# define your gauss points\n",
    "qpoints=gausspoint2d(ngp=2)\n",
    "qpoints.creategausspoint()\n",
    "# for surface integration\n",
    "qpoints1d=gausspoint1d(ngp=2)\n",
    "qpoints1d.creategausspoint()\n",
    "\n",
    "#################################################\n",
    "### defines your physical parameters here\n",
    "#################################################\n",
    "D0=1.0e0\n",
    "j0=-0.02\n",
    "def D(x):\n",
    "    return D0*(1.0+0.8*np.sin(x*np.pi))\n",
    "def dD(x):\n",
    "    return D0*0.8*np.cos(x*np.pi)*np.pi\n",
    "    # return 2*x\n",
    "### for f\n",
    "def f(x):\n",
    "    return 1.0+10.0*x*x\n",
    "def df(x):\n",
    "    return 20.0*x\n",
    "\n",
    "# define your K and rhs, and other vectors\n",
    "nDofs=mymesh.nodes*1\n",
    "K   =np.zeros((nDofs,nDofs))\n",
    "R   =np.zeros(nDofs)\n",
    "C   =np.zeros(nDofs) # solution of current step\n",
    "Cold=np.zeros(nDofs) # solution of previous step\n",
    "penalty=1.0e16 # for dirichlet bc\n",
    "\n",
    "#################################################\n",
    "### for time stepping and NR iteration\n",
    "#################################################\n",
    "### for time stepping\n",
    "dt=2.0e-3     # time step size\n",
    "totalstep=50 # total time step\n",
    "### for NR iteration\n",
    "tolerance=1.0e-10;reltol=1.0e-15 # absolute tolerenace and relative tolerance\n",
    "maxiters=50\n",
    "rnorm=1.0;rnorm0=1.0;iters=0 # init\n",
    "converged=False\n",
    "\n",
    "### for result output\n",
    "result=ResultIO('nldiff')\n",
    "result.save2vtu(mesh=mymesh,solution=C,varnamelist=['c'],step=0) # for initial configuration\n",
    "interval=5 # output per 5-steps\n",
    "\n",
    "gradc=np.zeros(2)\n",
    "for step in range(1,totalstep):\n",
    "    C[:]=Cold[:]\n",
    "    iters=0;converged=False\n",
    "    # now we can do the newton-raphson iteration\n",
    "    print('step=%6d, dt=%14.5e, total time=%14.5e'%(step,dt,dt*step))\n",
    "    while iters<maxiters and (not converged):\n",
    "        K[:,:]=0.0;R[:]=0.0 # init\n",
    "        # do the element loop\n",
    "        for e in range(mymesh.elements):\n",
    "            elconn=mymesh.elementconn[e,:]\n",
    "            nodes=mymesh.nodecoords[elconn]\n",
    "            for gp in range(qpoints.ngp2):\n",
    "                xi =qpoints.gpcoords[gp,1]\n",
    "                eta=qpoints.gpcoords[gp,2]\n",
    "                w  =qpoints.gpcoords[gp,0]\n",
    "                shp_val,shp_grad,j=shp.calc(xi,eta,nodes[:,0],nodes[:,1]) # xi,eta,x,y\n",
    "                JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "                \n",
    "                ###########################################################\n",
    "                ### calc the physical quantities on current gauss point\n",
    "                ###########################################################\n",
    "                c   =0.0 # for the concentration on current gauss point\n",
    "                cold=0.0 # for the concentration of previous step\n",
    "                gradc[0]=0.0;gradc[1]=0.0 # gradient of concentration\n",
    "                x=0.0;y=0.0 # coordinates of current gauss point\n",
    "                for i in range(mymesh.nodesperelement):\n",
    "                    iInd=elconn[i] # global id\n",
    "                    c   +=shp_val[i]*C[iInd]\n",
    "                    cold+=shp_val[i]*Cold[iInd]\n",
    "                    gradc[0]+=shp_grad[i,0]*C[iInd]\n",
    "                    gradc[1]+=shp_grad[i,1]*C[iInd]\n",
    "                    x+=shp_val[i]*nodes[i,0]\n",
    "                    y+=shp_val[i]*nodes[i,1]\n",
    "                # here we use backward euler for the time integration\n",
    "                cdot=(c-cold)/(dt) # for the concentration 'rate'\n",
    "                dcdotdc=1.0/dt     # for \\frac{\\partial\\dot{c}}{\\partial c}\n",
    "                for i in range(mymesh.nodesperelement):\n",
    "                    iInd=elconn[i]\n",
    "                    # assemble to residual\n",
    "                    # R_{c}^{I}=\\int_{\\Omega}\\dot{c}N^{I}dV\n",
    "                    #          +\\int_{\\Omega}D(c)\\nabla c\\nabla N^{I}dV\n",
    "                    #          -\\int_{\\Omega}f(c)N^{I}dV\n",
    "                    #          +\\int_{\\partial\\Omega}j_{0}N^{I}dS\n",
    "                    R[iInd]+=cdot*shp_val[i]*JxW \\\n",
    "                            +D(c)*(gradc[0]*shp_grad[i,0]+gradc[1]*shp_grad[i,1])*JxW\\\n",
    "                            -f(c)*shp_val[i]*JxW\n",
    "                    # now we can calculate the K matrix\n",
    "                    for j in range(mymesh.nodesperelement):\n",
    "                        jInd=elconn[j]\n",
    "                        # K_{cc}^{IJ}&=-\\frac{\\partial R_{c}^{I}}{\\partial c^{J}}\\\\\n",
    "                        # =-\\int_{\\Omega}\\frac{\\partial\\dot{c}}{\\partial c}N^{J}N^{I}dV\\\\\n",
    "                        #  -\\int_{\\Omega}\\frac{\\partial D(c)}{\\partial c}N^{J}\\nabla c\\nabla N^{I} dV\n",
    "                        #  -\\int_{\\Omega}D(c)\\nabla N^{J}\\nabla N^{I}dV\\\\\n",
    "                        #  +\\int_{\\Omega}\\frac{\\partial f(c)}{\\partial c}N^{J}N^{I}dV\n",
    "                        K[iInd,jInd]+=(-1.0)*dcdotdc*shp_val[j]*shp_val[i]*JxW \\\n",
    "                                     +(-1.0)*dD(c)*shp_val[j]*(gradc[0]*shp_grad[i,0]+gradc[1]*shp_grad[i,1])*JxW \\\n",
    "                                     +(-1.0)*D(c)*(shp_grad[j,0]*shp_grad[i,0]+shp_grad[j,1]*shp_grad[i,1])*JxW \\\n",
    "                                     +(+1.0)*df(c)*shp_val[j]*shp_val[i]*JxW\n",
    "            # end-of-gauss-point-loop\n",
    "        # end-of-volume-element-loop\n",
    "        #################################################\n",
    "        ### now we need to apply the flux to 'bc' edge\n",
    "        #################################################\n",
    "#         ### here we apply the flux to 'right' edge\n",
    "#         bcelmt=mymesh.bcconn['right'] # we want to apply the flux on the right side\n",
    "#         for e in range(np.shape(bcelmt)[0]):\n",
    "#             elconn=bcelmt[e,:]\n",
    "#             nodes=mymesh.nodecoords[elconn,:]\n",
    "#             for gp in range(qpoints1d.ngp):\n",
    "#                 xi =qpoints1d.gpcoords[gp,1] # 1->xi, 2->eta\n",
    "#                 w  =qpoints1d.gpcoords[gp,0] # weight\n",
    "        \n",
    "#                 shp_val,shp_grad,j=shp1d.calc(xi,nodes[:,1]) # for left side, you need 'y'\n",
    "#                 JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "#                 # now we can calculate the surface flux\n",
    "#                 y=0.0\n",
    "#                 for i in range(len(elconn)):\n",
    "#                     y+=shp_val[i]*nodes[i,1]\n",
    "#                 for i in range(len(elconn)):\n",
    "#                     iInd=elconn[i]\n",
    "#                     # R_{c}^{I}=\\int_{\\Omega}\\dot{c}N^{I}dV\n",
    "#                     #          +\\int_{\\Omega}D(c)\\nabla c\\nabla N^{I}dV\n",
    "#                     #          -\\int_{\\Omega}f(c)N^{I}dV\n",
    "#                     #          +\\int_{\\partial\\Omega}j_{0}N^{I}dS\n",
    "#                     # assemble to global residual\n",
    "#                     # R[iInd]+=j0*shp_val[i]*JxW\n",
    "#                     R[iInd]+=j0*np.sin(y)*shp_val[i]*JxW\n",
    "        #         ### here we apply the flux to 'right' edge\n",
    "        bcelmt=mymesh.bcconn['top'] # we want to apply the flux on the right side\n",
    "        for e in range(np.shape(bcelmt)[0]):\n",
    "            elconn=bcelmt[e,:]\n",
    "            nodes=mymesh.nodecoords[elconn,:]\n",
    "            for gp in range(qpoints1d.ngp):\n",
    "                xi =qpoints1d.gpcoords[gp,1] # 1->xi, 2->eta\n",
    "                w  =qpoints1d.gpcoords[gp,0] # weight\n",
    "        \n",
    "                shp_val,shp_grad,j=shp1d.calc(xi,nodes[:,0]) # for left side, you need 'x'\n",
    "                JxW=j*w # the weight times the determinate of jacobian transformation\n",
    "                # now we can calculate the surface flux\n",
    "                x=0.0\n",
    "                for i in range(len(elconn)):\n",
    "                    x+=shp_val[i]*nodes[i,0]\n",
    "                for i in range(len(elconn)):\n",
    "                    iInd=elconn[i]\n",
    "                    # R_{c}^{I}=\\int_{\\Omega}\\dot{c}N^{I}dV\n",
    "                    #          +\\int_{\\Omega}D(c)\\nabla c\\nabla N^{I}dV\n",
    "                    #          -\\int_{\\Omega}f(c)N^{I}dV\n",
    "                    #          +\\int_{\\partial\\Omega}j_{0}N^{I}dS\n",
    "                    # assemble to global residual\n",
    "                    # R[iInd]+=j0*shp_val[i]*JxW\n",
    "                    R[iInd]+=j0*np.sin(x*np.pi*2.0)*shp_val[i]*JxW\n",
    "        # end-of-bc-element-loop\n",
    "        ###################################################\n",
    "        ### now your K and R are ready, we can solve Ax=b\n",
    "        ###################################################\n",
    "        dC=np.linalg.solve(K,R)\n",
    "        # update iteration info\n",
    "        iters+=1\n",
    "        C[:]+=dC[:]\n",
    "        rnorm=np.linalg.norm(R)\n",
    "        dunorm=np.linalg.norm(dC)\n",
    "        enorm=rnorm*dunorm\n",
    "        if iters==1:\n",
    "            rnorm0=rnorm\n",
    "            dunorm0=dunorm\n",
    "            enorm0=enorm\n",
    "        print('    iters=%3d, |R0 |=%14.5e, |R |=%14.5e'%(iters,rnorm0,rnorm))\n",
    "        # print('           |dU0|=%14.5e, |dU|=%14.5e'%(dunorm0,dunorm))\n",
    "        # print('           |E0 |=%14.5e, |E |=%14.5e'%(enorm0,enorm))\n",
    "        ##########################################\n",
    "        # convergence check\n",
    "        ##########################################\n",
    "        if rnorm<tolerance:\n",
    "            converged=True\n",
    "            break\n",
    "    ### end-of-NR-loop\n",
    "    ##################################\n",
    "    ### update solution info and output\n",
    "    ##################################\n",
    "    Cold[:]=C[:]\n",
    "    if step%interval==0 or step==totalstep-1:\n",
    "        result.save2vtu(mesh=mymesh,solution=C,varnamelist=['c'],step=step)\n",
    "        print('---------------------------------------------------------')\n",
    "### end-of-time-step-loop\n",
    "print('simulation is done!')"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}