diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 68541b5..7d0efab 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -1,4 +1,4 @@
-name: Test notebook execution
+name: Test and deploy notebooks
on:
pull_request:
@@ -6,7 +6,7 @@ on:
jobs:
build:
- name: Test notebook execution
+ name: Test and deploy notebooks
runs-on: self-hosted
container:
image: firedrakeproject/firedrake:latest
@@ -19,7 +19,7 @@ jobs:
- uses: actions/checkout@v3
- name: Load pyrol wheel from cache
- uses: actions/cache@v3
+ uses: actions/cache/restore@v3
id: cache-pyrol
with:
path: /tmp/wheels
@@ -28,16 +28,47 @@ jobs:
- name: Build pyrol wheel
if: steps.cache-pyrol.outputs.cache-hit != 'true'
run: |
- /home/firedrake/firedrake/bin/pip wheel -r git+https://github.com/angus-g/pyrol@rol-2.0-checkpointing -w /tmp/wheels
+ /home/firedrake/firedrake/bin/pip wheel git+https://github.com/angus-g/pyrol@rol-2.0-checkpointing -w /tmp/wheels
- - name: Install Python dependencies
+ - name: Save pyrol to cache
+ if: steps.cache-pyrol.outputs.cache-hit != 'true'
+ uses: actions/cache/save@v3
+ with:
+ path: /tmp/wheels
+ key: pyrol
+
+ - name: Install dependencies
run: |
+ sudo apt update
+ sudo apt install -y libgl1-mesa-glx xvfb
. /home/firedrake/firedrake/bin/activate
python3 -m pip install --no-index -f /tmp/wheels pyroltrilinos
- python3 -m pip install nbval
- python3 -m pip install git+https://github.com/g-adopt/g-adopt@adjoint-cases-testing
+ python3 -m pip install nbval pyvista
+ python3 -m pip install git+https://github.com/g-adopt/g-adopt
- name: Run test
run: |
. /home/firedrake/firedrake/bin/activate
+ export DISPLAY=:99
+ export PYVISTA_OFF_SCREEN=true
+ Xvfb $DISPLAY -screen 0 1024x768x24 > /dev/null 2>&1 &
+ sleep 3
pytest --nbval *.ipynb
+
+ - name: Build Jupyter Book
+ run: |
+ . /home/firedrake/firedrake/bin/activate
+ python3 -m pip install jupyter-book
+ for f in *.ipynb; do ln -s ../$f docs; done
+ ln -s ../stokes-control.msh docs
+ export DISPLAY=:99
+ export PYVISTA_OFF_SCREEN=true
+ Xvfb $DISPLAY -screen 0 1024x768x24 > /dev/null 2>&1 &
+ sleep 3
+ jupyter-book build docs
+
+ - name: Archive built book
+ uses: actions/upload-artifact@v3
+ with:
+ name: tutorials-html
+ path: docs/_build/html
diff --git a/01-spd-helmholtz.ipynb b/01-spd-helmholtz.ipynb
index 5134ffa..dbef946 100644
--- a/01-spd-helmholtz.ipynb
+++ b/01-spd-helmholtz.ipynb
@@ -7,14 +7,14 @@
"outputs": [],
"source": [
"# This magic makes plots appear in the browser\n",
- "%matplotlib notebook\n",
+ "%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# Load Firedrake on Colab\n",
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
" import firedrake"
]
},
@@ -52,7 +52,9 @@
{
"cell_type": "markdown",
"metadata": {
- "collapsed": true
+ "jupyter": {
+ "outputs_hidden": true
+ }
},
"source": [
"$$ \\int_\\Omega \\nabla u\\cdot\\nabla v + uv\\ \\mathrm{d}x = \\int_\\Omega\n",
@@ -129,7 +131,7 @@
"metadata": {},
"outputs": [],
"source": [
- "from firedrake import *"
+ "from gadopt import *"
]
},
{
@@ -187,10 +189,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"# We test the output of the notebooks with out continuous integration system to\n",
"# make sure they run.\n",
"# Unfortunately we can't compare the plots from run to run, so the above line\n",
@@ -312,10 +317,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"fig, axes = plt.subplots()\n",
"collection = tripcolor(uh, axes=axes, cmap='coolwarm')\n",
"fig.colorbar(collection);"
@@ -340,10 +348,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"difference = assemble(interpolate(u_exact, V) - uh)\n",
"fig, axes = plt.subplots()\n",
"collection = tripcolor(difference, axes=axes, cmap='coolwarm')\n",
@@ -354,17 +365,17 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Exercises\n",
+ "## Exercises\n",
"\n",
- "## Exercise 1: \n",
+ "### Exercise 1: \n",
"\n",
- "### 1a: use a higher approximation degree\n",
+ "#### 1a: use a higher approximation degree\n",
"\n",
"Solve the same problem, only this time, use a piecewise quadratic approximation space.\n",
"\n",
"- Hint: check the help for `FunctionSpace` to see how to specify the degree.\n",
"\n",
- "### 1b: use a quadrilateral mesh\n",
+ "#### 1b: use a quadrilateral mesh\n",
"\n",
"Solve the same problem, but using quadrilateral, rather than triangular, cells.\n",
"\n",
@@ -383,14 +394,14 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Exercise 2: convergence of the method\n",
+ "### Exercise 2: convergence of the method\n",
"For solutions with sufficient smoothness (like the choice we have here), this method with a piecewise linear approximation space should converge in the $L_2$ error with rate $\\mathcal{O}(h^{-2})$, where $h$ is the typical mesh spacing. Confirm this for the example in question by computing the $L_2$ error in the solution for a sequence of finer and finer meshes.\n",
"\n",
"- Hint 1: You can compute errors using [errornorm](http://firedrakeproject.org/firedrake.html#firedrake.norms.errornorm)\n",
"- Hint 2: If the error is $\\mathcal{O}(h^{-2})$ then $\\log_2 (e_H/e_h) \\approx 2$.\n",
" The Python module `math` contains an implementation of `log`.\n",
" \n",
- "### What works better? Mesh refinement, or increasing the approximation degree?\n",
+ "#### What works better? Mesh refinement, or increasing the approximation degree?\n",
"\n",
"Instead of (or as well as!) refining the mesh, we can increase the degree of the approximating polynomial space.\n",
"\n",
@@ -430,7 +441,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Part II: inhomogeneous Neuman conditions\n",
+ "## Part II: inhomogeneous Neuman conditions\n",
"\n",
"Let us recall again the statement of our problem. We seek $u \\in V$ satisfying\n",
"\n",
@@ -449,7 +460,9 @@
{
"cell_type": "markdown",
"metadata": {
- "collapsed": true
+ "jupyter": {
+ "outputs_hidden": true
+ }
},
"source": [
"As previously, we introduce a test function $v \\in V$, multiply and integrate. After integrating by parts, we obtain\n",
@@ -548,11 +561,6 @@
}
],
"metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
"language_info": {
"codemirror_mode": {
"name": "ipython",
@@ -563,9 +571,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.10.2"
+ "version": "3.11.4"
}
},
"nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 4
}
diff --git a/02-poisson.ipynb b/02-poisson.ipynb
index 203ade1..68278c4 100644
--- a/02-poisson.ipynb
+++ b/02-poisson.ipynb
@@ -10,7 +10,9 @@
{
"cell_type": "markdown",
"metadata": {
- "collapsed": true
+ "jupyter": {
+ "outputs_hidden": true
+ }
},
"source": [
"Let's move on from the Helmholtz problem to Poisson:\n",
@@ -49,15 +51,15 @@
"metadata": {},
"outputs": [],
"source": [
- "%matplotlib notebook\n",
+ "%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"# Load Firedrake on Colab\n",
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
" import firedrake\n",
- "from firedrake import *\n",
+ "from gadopt import *\n",
"mesh = UnitSquareMesh(10, 10)\n",
"V = FunctionSpace(mesh, \"Lagrange\", 1)"
]
@@ -159,10 +161,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"fig, axes = plt.subplots()\n",
"collection = tripcolor(uh, axes=axes)\n",
"fig.colorbar(collection);"
@@ -172,7 +177,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Exercises\n",
+ "## Exercises\n",
"\n",
"Most of the time, we don't want to impose the same Dirichlet condition everywhere. Instead of solving with homogeneous Dirichlet conditions everywhere, solve the following problem.\n",
"\n",
@@ -199,11 +204,6 @@
}
],
"metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
"language_info": {
"codemirror_mode": {
"name": "ipython",
@@ -214,9 +214,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.11.4"
}
},
"nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 4
}
diff --git a/03-elasticity.ipynb b/03-elasticity.ipynb
index 1c38933..936b697 100644
--- a/03-elasticity.ipynb
+++ b/03-elasticity.ipynb
@@ -47,15 +47,15 @@
"metadata": {},
"outputs": [],
"source": [
- "%matplotlib notebook\n",
+ "%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"# Load Firedrake on Colab\n",
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
" import firedrake\n",
- "from firedrake import *\n",
+ "from gadopt import *\n",
"length = 1\n",
"width = 0.2\n",
"mesh = RectangleMesh(40, 20, length, width)"
@@ -185,10 +185,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"fig, axes = plt.subplots()\n",
"triplot(displaced_mesh, axes=axes)\n",
"axes.set_aspect(\"equal\");"
@@ -198,7 +201,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Exercises\n",
+ "## Exercises\n",
"\n",
"Modify the problem so that the material density $\\rho$ is not constant, but rather varies in space. For example, you could try setting $\\rho(x, y) = 0.01 + xy$. That is, the material gets denser the futher away from the clamped end you get.\n",
"\n",
@@ -217,7 +220,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Solving bigger problems\n",
+ "## Solving bigger problems\n",
"\n",
"Up to now, we've only really solved quite small problems, and therefore haven't really had to worry about tuning the solver. As we increase the size of the problem we're solving, the direct solver approach will no longer be good enough. Firedrake uses [PETSc](http://www.mcs.anl.gov/petsc) to provide solvers, and uses PETSc solver parameters to control them.\n",
"\n",
@@ -240,12 +243,7 @@
" f = as_vector([0, -rho*g])\n",
" mu = Constant(1)\n",
" lambda_ = Constant(0.25)\n",
- " Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor\n",
- " def epsilon(u):\n",
- " return 0.5*(grad(u) + grad(u).T)\n",
- "\n",
- " def sigma(u):\n",
- " return lambda_*div(u)*Id + 2*mu*epsilon(u) \n",
+ " Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor \n",
" \n",
" bc = DirichletBC(V, Constant([0, 0]), 1)\n",
" u = TrialFunction(V)\n",
@@ -269,10 +267,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "raises-exception"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_RAISES_EXCEPTION\n",
"uh = solve_elasticity(100, 100, options={\"ksp_max_it\": 100, \"ksp_type\": \"cg\", \"pc_type\": \"none\"})"
]
},
@@ -293,6 +294,7 @@
" \"ksp_max_it\": 100, \n",
" \"pc_type\": \"gamg\",\n",
" \"mg_levels_pc_type\": \"sor\",\n",
+ " \"mat_type\": \"aij\",\n",
" \"ksp_converged_reason\": None})"
]
},
@@ -320,11 +322,7 @@
" mu = Constant(1)\n",
" lambda_ = Constant(0.25)\n",
" Id = Identity(mesh.geometric_dimension()) # 2x2 Identity tensor\n",
- " def epsilon(u):\n",
- " return 0.5*(grad(u) + grad(u).T)\n",
- "\n",
- " def sigma(u):\n",
- " return lambda_*div(u)*Id + 2*mu*epsilon(u) \n",
+ " \n",
" bc = DirichletBC(V, Constant([0, 0]), 1)\n",
" u = TrialFunction(V)\n",
" v = TestFunction(V)\n",
@@ -351,7 +349,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "With this done, the problem is solved in a reasonably small number of Krylov iterations."
+ "With this done, the problem is solved in a reasonably small number of Krylov iterations. This does require that the GAMG smoothing is switched back to the default of Jacobi, rather than SOR as above."
]
},
{
@@ -363,7 +361,7 @@
"uh = solve_elasticity(200, 200, options={\"ksp_type\": \"cg\", \n",
" \"ksp_max_it\": 100, \n",
" \"pc_type\": \"gamg\",\n",
- " \"mg_levels_pc_type\": \"sor\",\n",
+ " \"mat_type\": \"aij\",\n",
" \"ksp_converged_reason\": None})"
]
},
@@ -385,11 +383,6 @@
}
],
"metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
"language_info": {
"codemirror_mode": {
"name": "ipython",
@@ -400,9 +393,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.11.4"
}
},
"nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 4
}
diff --git a/04-burgers.ipynb b/04-burgers.ipynb
index bfccfb8..ade8e18 100644
--- a/04-burgers.ipynb
+++ b/04-burgers.ipynb
@@ -55,9 +55,9 @@
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
" import firedrake\n",
- "from firedrake import *\n",
+ "from gadopt import *\n",
"\n",
"n = 100\n",
"mesh = PeriodicIntervalMesh(n, length=2)\n",
@@ -196,10 +196,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"import matplotlib.pyplot as plt\n",
"fig, axes = plt.subplots()\n",
"axes.set_ylim((-1., 1.))\n",
@@ -218,10 +221,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"from matplotlib.animation import FuncAnimation\n",
"\n",
"def animate(u):\n",
@@ -344,11 +350,6 @@
}
],
"metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
"language_info": {
"codemirror_mode": {
"name": "ipython",
@@ -359,258 +360,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
- },
- "widgets": {
- "application/vnd.jupyter.widget-state+json": {
- "state": {
- "483f8c5b55a94cb5b9b3c6223be10e49": {
- "model_module": "@jupyter-widgets/controls",
- "model_module_version": "1.5.0",
- "model_name": "VBoxModel",
- "state": {
- "_dom_classes": [
- "widget-interact"
- ],
- "_model_module": "@jupyter-widgets/controls",
- "_model_module_version": "1.5.0",
- "_model_name": "VBoxModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/controls",
- "_view_module_version": "1.5.0",
- "_view_name": "VBoxView",
- "box_style": "",
- "children": [
- "IPY_MODEL_8f43ada5d62046f7a6a50fbcc6587acf",
- "IPY_MODEL_ff45eb9314f84c0b99d6cbc4dafaeea7"
- ],
- "layout": "IPY_MODEL_c987f4ac87ac4cae9b51332042a41dd2"
- }
- },
- "8f43ada5d62046f7a6a50fbcc6587acf": {
- "model_module": "@jupyter-widgets/controls",
- "model_module_version": "1.5.0",
- "model_name": "IntSliderModel",
- "state": {
- "_dom_classes": [],
- "_model_module": "@jupyter-widgets/controls",
- "_model_module_version": "1.5.0",
- "_model_name": "IntSliderModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/controls",
- "_view_module_version": "1.5.0",
- "_view_name": "IntSliderView",
- "continuous_update": true,
- "description": "index",
- "description_tooltip": null,
- "disabled": false,
- "layout": "IPY_MODEL_e34369af8020497289732bb687d2e8e2",
- "max": 50,
- "min": 0,
- "orientation": "horizontal",
- "readout": true,
- "readout_format": "d",
- "step": 1,
- "style": "IPY_MODEL_eb5e810f3cf7471c8f32d7e6ef99b66b",
- "value": 0
- }
- },
- "c987f4ac87ac4cae9b51332042a41dd2": {
- "model_module": "@jupyter-widgets/base",
- "model_module_version": "1.2.0",
- "model_name": "LayoutModel",
- "state": {
- "_model_module": "@jupyter-widgets/base",
- "_model_module_version": "1.2.0",
- "_model_name": "LayoutModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/base",
- "_view_module_version": "1.2.0",
- "_view_name": "LayoutView",
- "align_content": null,
- "align_items": null,
- "align_self": null,
- "border": null,
- "bottom": null,
- "display": null,
- "flex": null,
- "flex_flow": null,
- "grid_area": null,
- "grid_auto_columns": null,
- "grid_auto_flow": null,
- "grid_auto_rows": null,
- "grid_column": null,
- "grid_gap": null,
- "grid_row": null,
- "grid_template_areas": null,
- "grid_template_columns": null,
- "grid_template_rows": null,
- "height": null,
- "justify_content": null,
- "justify_items": null,
- "left": null,
- "margin": null,
- "max_height": null,
- "max_width": null,
- "min_height": null,
- "min_width": null,
- "object_fit": null,
- "object_position": null,
- "order": null,
- "overflow": null,
- "overflow_x": null,
- "overflow_y": null,
- "padding": null,
- "right": null,
- "top": null,
- "visibility": null,
- "width": null
- }
- },
- "e34369af8020497289732bb687d2e8e2": {
- "model_module": "@jupyter-widgets/base",
- "model_module_version": "1.2.0",
- "model_name": "LayoutModel",
- "state": {
- "_model_module": "@jupyter-widgets/base",
- "_model_module_version": "1.2.0",
- "_model_name": "LayoutModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/base",
- "_view_module_version": "1.2.0",
- "_view_name": "LayoutView",
- "align_content": null,
- "align_items": null,
- "align_self": null,
- "border": null,
- "bottom": null,
- "display": null,
- "flex": null,
- "flex_flow": null,
- "grid_area": null,
- "grid_auto_columns": null,
- "grid_auto_flow": null,
- "grid_auto_rows": null,
- "grid_column": null,
- "grid_gap": null,
- "grid_row": null,
- "grid_template_areas": null,
- "grid_template_columns": null,
- "grid_template_rows": null,
- "height": null,
- "justify_content": null,
- "justify_items": null,
- "left": null,
- "margin": null,
- "max_height": null,
- "max_width": null,
- "min_height": null,
- "min_width": null,
- "object_fit": null,
- "object_position": null,
- "order": null,
- "overflow": null,
- "overflow_x": null,
- "overflow_y": null,
- "padding": null,
- "right": null,
- "top": null,
- "visibility": null,
- "width": null
- }
- },
- "eaddcef7c5f24f0bae9a44e7536afc39": {
- "model_module": "@jupyter-widgets/base",
- "model_module_version": "1.2.0",
- "model_name": "LayoutModel",
- "state": {
- "_model_module": "@jupyter-widgets/base",
- "_model_module_version": "1.2.0",
- "_model_name": "LayoutModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/base",
- "_view_module_version": "1.2.0",
- "_view_name": "LayoutView",
- "align_content": null,
- "align_items": null,
- "align_self": null,
- "border": null,
- "bottom": null,
- "display": null,
- "flex": null,
- "flex_flow": null,
- "grid_area": null,
- "grid_auto_columns": null,
- "grid_auto_flow": null,
- "grid_auto_rows": null,
- "grid_column": null,
- "grid_gap": null,
- "grid_row": null,
- "grid_template_areas": null,
- "grid_template_columns": null,
- "grid_template_rows": null,
- "height": null,
- "justify_content": null,
- "justify_items": null,
- "left": null,
- "margin": null,
- "max_height": null,
- "max_width": null,
- "min_height": null,
- "min_width": null,
- "object_fit": null,
- "object_position": null,
- "order": null,
- "overflow": null,
- "overflow_x": null,
- "overflow_y": null,
- "padding": null,
- "right": null,
- "top": null,
- "visibility": null,
- "width": null
- }
- },
- "eb5e810f3cf7471c8f32d7e6ef99b66b": {
- "model_module": "@jupyter-widgets/controls",
- "model_module_version": "1.5.0",
- "model_name": "SliderStyleModel",
- "state": {
- "_model_module": "@jupyter-widgets/controls",
- "_model_module_version": "1.5.0",
- "_model_name": "SliderStyleModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/base",
- "_view_module_version": "1.2.0",
- "_view_name": "StyleView",
- "description_width": "",
- "handle_color": null
- }
- },
- "ff45eb9314f84c0b99d6cbc4dafaeea7": {
- "model_module": "@jupyter-widgets/output",
- "model_module_version": "1.0.0",
- "model_name": "OutputModel",
- "state": {
- "_dom_classes": [],
- "_model_module": "@jupyter-widgets/output",
- "_model_module_version": "1.0.0",
- "_model_name": "OutputModel",
- "_view_count": null,
- "_view_module": "@jupyter-widgets/output",
- "_view_module_version": "1.0.0",
- "_view_name": "OutputView",
- "layout": "IPY_MODEL_eaddcef7c5f24f0bae9a44e7536afc39",
- "msg_id": "",
- "outputs": []
- }
- }
- },
- "version_major": 2,
- "version_minor": 0
- }
+ "version": "3.11.4"
}
},
"nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 4
}
diff --git a/05-pde-constrained-optimisation.ipynb b/05-pde-constrained-optimisation.ipynb
index 1417799..76370b3 100644
--- a/05-pde-constrained-optimisation.ipynb
+++ b/05-pde-constrained-optimisation.ipynb
@@ -9,17 +9,7 @@
"*This example is modified from the equivalent [`dolfin-adjoint` demo](http://www.dolfin-adjoint.org/en/latest/documentation/stokes-bc-control/stokes-bc-control.html)*\n",
"\n",
"\n",
- "In this example, we will look at how to use `firedrake-adjoint` to optimise for strong (Dirichlet) conditions in a steady problem. `firedrake-adjoint` is a thin compatibility layer for the [`dolfin-adjoint` package](http://www.dolfin-adjoint.org/en/latest/), a python package to **automatically derive the discrete adjoint and tangent linear models** of forward problems written using Firedrake.\n",
- "\n",
- "## Installing necessary dependencies\n",
- "\n",
- "For the minimisation, we will need `scipy`, which is installed with\n",
- "\n",
- "```sh\n",
- "pip install scipy\n",
- "```\n",
- "\n",
- "In the activated firedrake virtualenv."
+ "In this example, we will look at how to use `firedrake-adjoint` to optimise for strong (Dirichlet) conditions in a steady problem. `firedrake-adjoint` is a thin compatibility layer for the [`dolfin-adjoint` package](http://www.dolfin-adjoint.org/en/latest/), a python package to **automatically derive the discrete adjoint and tangent linear models** of forward problems written using Firedrake."
]
},
{
@@ -35,22 +25,21 @@
"metadata": {},
"outputs": [],
"source": [
- "%matplotlib notebook\n",
+ "%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"# Load Firedrake on Colab\n",
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
- " import firedrake\n",
- "from firedrake import *"
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " import firedrake"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "The next step is new, we now import the `firedrake_adjoint` package, which *overrides* much of the Firedrake interface so that an *annotated tape* for adjoint calculations can be built automatically."
+ "The next step is new, we now import the `gadopt.inverse` package, which *overrides* much of the Firedrake interface through `firedrake.adjoint` so that an *annotated tape* for adjoint calculations can be built automatically."
]
},
{
@@ -59,7 +48,8 @@
"metadata": {},
"outputs": [],
"source": [
- "from firedrake_adjoint import *"
+ "from gadopt import *\n",
+ "from gadopt.inverse import *"
]
},
{
@@ -101,6 +91,15 @@
"This problem setup requires a mesh that is more complex than the built in ones Firedrake provides. Instead, it was created with [Gmsh](http://gmsh.info). It is loaded by using the `Mesh` constructor, passing the filename of the mesh in question."
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "!wget https://raw.githubusercontent.com/g-adopt/tutorials/main/stokes-control.msh"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -120,10 +119,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"fig, axes = plt.subplots()\n",
"triplot(mesh, axes=axes)\n",
"axes.axis(\"off\")\n",
@@ -236,10 +238,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"u_init, p_init = w.split()\n",
"\n",
"fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)\n",
@@ -296,10 +301,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"fig, axes = plt.subplots()\n",
"arrows = quiver(g_opt, axes=axes, scale=3)\n",
"axes.set_aspect(\"equal\")\n",
@@ -347,10 +355,13 @@
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# NBVAL_IGNORE_OUTPUT\n",
"u_opt, p_opt = w_opt.split()\n",
"\n",
"fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)\n",
@@ -396,11 +407,6 @@
}
],
"metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
"language_info": {
"codemirror_mode": {
"name": "ipython",
@@ -411,9 +417,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.4"
+ "version": "3.11.4"
}
},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/06-GD-2D-convection.ipynb b/06-GD-2D-convection.ipynb
index 1e51158..3399f14 100644
--- a/06-GD-2D-convection.ipynb
+++ b/06-GD-2D-convection.ipynb
@@ -5,12 +5,11 @@
"metadata": {},
"source": [
"# Idealised 2-D mantle convection problem in a square box\n",
- "We next turn our attention to the slow creeping motion of Earth’s mantle \n",
+ "In this tutorial, we examine the slow creeping motion of Earth’s mantle\n",
"over geological timescales. We start with an idealised 2-D problem. \n",
"\n",
"## Governing equations\n",
- "The equations governing mantle convection\n",
- "are derived from the conservation laws of mass, momentum and energy.\n",
+ "The equations governing mantle convection are derived from the conservation laws of mass, momentum and energy.\n",
"The simplest mathematical formulation assumes incompressibility and\n",
"the Boussinesq approximation (e.g. McKenzie et al., 1973), under which the\n",
"non–dimensional momentum and continuity equations are given by: \n",
@@ -23,7 +22,7 @@
"\n",
"$$Ra0 = \\frac{\\rho_0 \\alpha \\Delta T g d^3}{\\mu_0 \\kappa}$$\n",
"\n",
- "Here, $\\rho_0$ denotes reference density, $\\alpha$ the thermal expansion \n",
+ "Here, $\\rho_0$ denotes reference density, $\\alpha$ the thermal expansion\n",
"coefficient, $\\Delta T$ the characteristic temperature change across the\n",
"domain, $g$ the gravitational acceleration, $d$ the characteristic length,\n",
"$\\mu_0$ the reference dynamic viscosity and $\\kappa$ the thermal\n",
@@ -33,17 +32,21 @@
" When simulating incompressible flow, the full stress tensor, $\\sigma$, is decomposed into deviatoric and volumetric components:\n",
"$$ \\sigma = \\tau - p I,$$\n",
"where $\\tau$ is the deviatoric stress tensor, $p$ is dynamic pressure and $I$ is the identity matrix. Substituting this into the first equation presented above and utilizing the constitutive relation,\n",
+ "\n",
"$$\\tau = 2\\mu \\dot\\epsilon =\n",
" 2\\mu \\operatorname{sym}(\\nabla \\vec{u}) =\n",
" \\mu\\left[ \\nabla \\vec{u} + \\left(\\nabla \\vec{u}\\right)^T\\right] $$\n",
+ "\n",
"which relates the deviatoric stress tensor, $\\tau$, to the strain-rate tensor, $\\dot\\epsilon=\\operatorname{sym}(\\nabla \\vec{u})$, yields:\n",
+ "\n",
"$$ \\nabla \\cdot \\mu \\left[{\\nabla\\vec{u}} + \\left(\\nabla\\vec{u}\\right)^T\\right]\n",
" - \\nabla p + Ra_{0} T\\hat{\\vec{k}} = 0. $$\n",
- "The viscous flow problem can thus be posed in terms of pressure, $p$, velocity, $\\vec{u}$, and temperature, $T$. \n",
+ "\n",
+ "The viscous flow problem can thus be posed in terms of pressure, $p$, velocity, $\\vec{u}$, and temperature, $T$.\n",
"\n",
"The evolution of the thermal field is controlled by an advection-diffusion equation, where, for simplicity, we ignore internal heating:\n",
"$$ \\frac{\\partial T}{\\partial t} + \\vec{u}\\cdot \\nabla T - \\nabla \\cdot \\left(\\kappa \\nabla T\\right) = 0 $$\n",
- "These governing equations are sufficient to solve for the three unknowns, together with adequate boundary and initial conditions. \n",
+ "These governing equations are sufficient to solve for the three unknowns, together with adequate boundary and initial conditions.\n",
"\n",
"## Weak formulation \n",
"\n",
@@ -52,14 +55,14 @@
"temperature T , and also contain the test functions v, w and q. The weak form is then obtained by multiplying these equations\n",
"with the test functions and integrating over the domain $\\Omega$, \n",
"\n",
- "$$\\int_\\Omega (\\nabla \\vec{v})\\colon \\mu \\left[ \\nabla \\vec{u} + \\left( \\nabla \\vec{u} \\right)^T\\right] \\ dx \n",
+ "$$\\int_\\Omega (\\nabla \\vec{v})\\colon \\mu \\left[ \\nabla \\vec{u} + \\left( \\nabla \\vec{u} \\right)^T\\right] \\ dx\n",
" - \\int_{\\Omega} \\left( \\nabla \\cdot \\vec{v}\\right)\\ p \\ dx\n",
" - \\int_{\\Omega} Ra_0\\ T\\ \\vec{v}\\cdot\\hat{k} \\ dx = 0 \\ \\text{ for all } v\\in V,$$\n",
"\n",
"$$ \\int_\\Omega w \\nabla \\cdot \\vec{u} \\ dx\\ \\text{ for all } v\\in V,$$\n",
"\n",
"$$ \\int_\\Omega q\\frac{\\partial T}{\\partial t} \\ dx\n",
- " + \\int_\\Omega q \\vec{u}\\cdot \\nabla T \\ dx \n",
+ " + \\int_\\Omega q \\vec{u}\\cdot \\nabla T \\ dx\n",
" + \\int_\\Omega \\left(\\nabla q\\right) \\cdot \\left(\\kappa \\nabla T\\right) \\ dx = 0 \\text{ for all } q\\in Q.$$\n",
"\n",
"Note that we have integrated by parts the viscosity and pressure gradient terms in the Stokes equations, and the diffusion term in the energy equation, but have omitted the corresponding boundary terms.\n",
@@ -67,25 +70,29 @@
"## Solution procedure \n",
"\n",
"For temporal integration, we apply a simple $\\theta$ scheme to the energy equation:\n",
+ "\n",
"$$\n",
" F_{\\text{energy}}(q; T^{n+1}) :=\n",
" \\int_\\Omega q \\frac{T^{n+1} - T^n}{\\Delta t} dx\n",
" + \\int_\\Omega q\\vec{u}\\cdot\\nabla T^{n+\\theta} dx\n",
- " + \\int_\\Omega \\left(\\nabla q\\right)\\cdot \\left(\\kappa \\nabla \n",
- " T^{n+\\theta}\\right) dx = 0 \n",
+ " + \\int_\\Omega \\left(\\nabla q\\right)\\cdot \\left(\\kappa \\nabla\n",
+ " T^{n+\\theta}\\right) dx = 0\n",
" \\text{ for all } q\\in Q,\n",
"$$\n",
+ "\n",
"where\n",
"$$\n",
" T^{n+\\theta} = \\theta T^{n+1} + (1-\\theta) T^n\n",
"$$\n",
+ "\n",
"is interpolated between the temperature solutions $T^n$ and $T^{n+1}$ at the\n",
- "beginning and end of the $n+1$-th time step using a parameter $0\\leq\\theta\\leq 1$. \n",
- "In this example we use a Crank-Nicolson scheme, where $\\theta = 0.5$. \n",
+ "beginning and end of the $n+1$-th time step using a parameter $0\\leq\\theta\\leq 1$.\n",
+ "In this example we use a Crank-Nicolson scheme, where $\\theta = 0.5$.\n",
"\n",
"To simplify we will solve for velocity and pressure, $\\vec{u}$ and $p$, in a separate step before solving for the new temperature $T^{n+1}$. Since these weak equations need to hold for all test\n",
"functions $\\vec{v}\\in V$ and $w\\in W$ we can equivalently write, using a single\n",
"residual functional $F_{\\text{Stokes}}$:\n",
+ "\n",
"$$\n",
" F_{\\text{Stokes}}(\\vec{v},w; \\vec{u}, p) =\n",
" \\int_\\Omega \\left(\\nabla \\vec{v}\\right) \\colon \\mu \\left[{\\nabla\\vec{u}}\n",
@@ -94,8 +101,9 @@
" - \\int_\\Omega Ra_{0} T\\vec{v}\\cdot\\hat{\\vec{k}} dx\n",
" - \\int_\\Omega w \\nabla \\cdot \\vec{u} dx = 0\n",
" \\text{ for all } \\vec{v}\\in V,\n",
- " w\\in W, \n",
+ " w\\in W,\n",
"$$\n",
+ "\n",
"where we have multiplied the continuity equation with $-1$ to ensure symmetry\n",
"between the $\\nabla p$ and $\\nabla\\cdot u$ terms. This combined weak form that we simultaneously solve for a velocity $u\\in V$ and pressure $p\\in W$ is referred to as a mixed problem, and the combined solution $(u, p)$ is said to be found in the mixed function space $V\\oplus W$.\n",
"\n",
@@ -106,15 +114,6 @@
"## Implementation"
]
},
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "!sudo mv /etc/apt/sources.list.d/cuda.list{,.disabled}"
- ]
- },
{
"cell_type": "code",
"execution_count": null,
@@ -122,17 +121,29 @@
"outputs": [],
"source": [
"# This magic makes plots appear in the browser\n",
- "%matplotlib notebook\n",
+ "%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
+ "import os\n",
+ "os.environ[\"GADOPT_LOGLEVEL\"] = \"WARNING\"\n",
+ "\n",
"# Load Firedrake on Colab\n",
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
" import firedrake"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from gadopt import *"
+ ]
+ },
{
"cell_type": "markdown",
"metadata": {},
@@ -141,12 +152,12 @@
"for velocity and pressure, with Q2 elements for temperature. Firedrake user code\n",
"is written in Python, so the first step, illustrated above, is to import the Firedrake module.\n",
"\n",
- "We next need a mesh: for simple domains such as the unit square, \n",
- "Firedrake provides built-in meshing functions. As such, the following code \n",
- "defines the mesh, with 20 quadrilateral elements in x and y directions. We also tag boundary ID's. \n",
+ "We next need a mesh: for simple domains such as the unit square,\n",
+ "Firedrake provides built-in meshing functions. As such, the following code\n",
+ "defines the mesh, with 20 quadrilateral elements in x and y directions. We also tag boundary IDs.\n",
"Boundaries are automatically tagged by the built-in meshes supported by Firedrake. For the `UnitSquareMesh`\n",
"being used here, tag 1 corresponds to the plane $x=0$; 2 to $x=1$; 3 to $y=0$; and 4 to\n",
- "$y=1$. We name these left, right, bottom and top, respectively. "
+ "$y=1$. We name these `left`, `right`, `bottom` and `top`, respectively."
]
},
{
@@ -155,7 +166,6 @@
"metadata": {},
"outputs": [],
"source": [
- "from firedrake import *\n",
"mesh = UnitSquareMesh(20, 20, quadrilateral=True)\n",
"left, right, bottom, top = 1, 2, 3, 4 # Boundary IDs"
]
@@ -164,9 +174,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "We also need function spaces, which is achieved by associating the mesh with \n",
+ "We also need function spaces, which is achieved by associating the mesh with\n",
"the relevant finite element: V , W and Q are symbolic variables\n",
- "representing function spaces. They also contain the function space’s \n",
+ "representing function spaces. They also contain the function space’s\n",
"computational implementation, recording the association of degrees of freedom\n",
"with the mesh and pointing to the finite element basis."
]
@@ -177,9 +187,9 @@
"metadata": {},
"outputs": [],
"source": [
- "V = VectorFunctionSpace(mesh, family=\"CG\", degree=2) # Velocity function space (vector)\n",
- "W = FunctionSpace(mesh, family=\"CG\", degree=1) # Pressure function space (scalar)\n",
- "Q = FunctionSpace(mesh, family=\"CG\", degree=2) # Temperature function space (scalar)"
+ "V = VectorFunctionSpace(mesh, \"CG\", degree=2) # Velocity function space (vector)\n",
+ "W = FunctionSpace(mesh, \"CG\", degree=1) # Pressure function space (scalar)\n",
+ "Q = FunctionSpace(mesh, \"CG\", degree=2) # Temperature function space (scalar)"
]
},
{
@@ -187,7 +197,7 @@
"metadata": {},
"source": [
"Function spaces can be combined in the natural way to create mixed function spaces,\n",
- "combining the velocity and pressure function spaces to form a function space for the mixed Stokes problem, Z."
+ "combining the velocity and pressure spaces to form a function space for the mixed Stokes problem, Z."
]
},
{
@@ -203,8 +213,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "Test functions, v, w and q, are subsequently defined and we also specify functions to hold our solutions: z in the mixed function space, noting that a symbolic representation of the two parts – velocity and pressure – is obtained with\n",
- "split, and Told and Tnew, required for the Crank-Nicolson scheme used for temporal discretisation in our energy equation, where $T_\\theta$ is also defined."
+ "Test functions, v, w and q, are subsequently defined and we also specify functions to hold our solutions: z in the mixed function space, noting that a symbolic representation of the two parts – velocity and pressure – is obtained with `split`."
]
},
{
@@ -217,23 +226,21 @@
"q = TestFunction(Q)\n",
"z = Function(Z) # a field over the mixed function space Z.\n",
"u, p = split(z) # returns symbolic UFL expression for u and p\n",
- "Told, Tnew = Function(Q, name=\"OldTemp\"), Function(Q, name=\"NewTemp\")\n",
- "Ttheta = 0.5 * Tnew + 0.5 * Told # Temporal discretisation through Crank-Nicolson"
+ "T = Function(Q, name=\"Temperature\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Mantle convection is an initial and boundary-value problem. We assume the initial temperature distribution to be prescribed by \n",
+ "Mantle convection is an initial and boundary-value problem. We assume the initial temperature distribution to be prescribed by \n",
"\n",
"$T(x,y) = (1-y) + 0.05\\ cos(\\pi x)\\ sin(\\pi y)$ \n",
"\n",
- "In the following code, we first obtain symbolic expressions for coordinates in the physical mesh and subsequently use these to initialize the old temperature field. \n",
+ "In the following code, we first obtain symbolic expressions for coordinates in the physical mesh and subsequently use these to initialize the old temperature field.\n",
"This is where Firedrake transforms a symbolic operation into a numerical computation for the first time: the\n",
"`interpolate` method generates C code that evaluates this expression at the nodes\n",
- "of $T_{\\text{old}}$, and immediately executes it to populate the values of $T_{\\text{old}}$. We\n",
- "initialize $T_{\\text{new}}$ with the values of $T_{\\text{old}}$."
+ "of $T$."
]
},
{
@@ -243,15 +250,14 @@
"outputs": [],
"source": [
"X = SpatialCoordinate(mesh)\n",
- "Told.interpolate(1.0 - X[1] + 0.05 * cos(pi * X[0]) * sin(pi * X[1]))\n",
- "Tnew.assign(Told)"
+ "T.interpolate(1.0 - X[1] + 0.05 * cos(pi * X[0]) * sin(pi * X[1]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We can visualise the initial temperature field using Firedrake's built-in plotting functionality. "
+ "We can visualise the initial temperature field using Firedrake's built-in plotting functionality."
]
},
{
@@ -261,41 +267,18 @@
"outputs": [],
"source": [
"fig, axes = plt.subplots()\n",
- "collection = tripcolor(Told, axes=axes, cmap='coolwarm')\n",
+ "collection = tripcolor(T, axes=axes, cmap='coolwarm')\n",
"fig.colorbar(collection);"
]
},
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": []
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "Important constants in this problem (Rayleigh Number, Ra; viscosity, μ; thermal diffusivity, κ), in addition to the\n",
- "constant timestep $\\Delta t$ and unit vector $\\hat{k}$, are next defined. We note that viscosity could also be a Function, if we\n",
- "wanted spatial variation, and will return to this in part 2 of this notebook below."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "Ra, mu, kappa, delta_t = Constant(1e4), Constant(1.0), Constant(1.0), Constant(1e-3)\n",
- "k = Constant((0, 1)) "
- ]
- },
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We are now in a position to define the variational problems expressed in our weak formulation above.\n",
- "We maintain the more general nonlinear residual form $F_{Stokes}(v, u) = 0$ and $F_{energy}(q, T)$ = 0, to allow\n",
- "for straightforward extension to nonlinear problems. Here we provide the symbolic expressions for $F_{Stokes}$ and $F_{Energy}$ in UFL: the resemblance to the mathematical formulation is immediately apparent. Integration over the domain is indicated by multiplication with `dx`."
+ "The Rayleigh number for this problem is defined in addition to the\n",
+ "initial timestep $\\Delta t$. The viscosity and thermal diffusivity are left at their default values (both = 1). We note that viscosity could also be a Function, if we wanted spatial variation, and will return to this in Part 2 of this notebook below.\n",
+ "\n",
+ "These constants are used to create an *Approximation* representing the physical setup of the problem (options include Boussinesq, Extended Boussinesq, Truncated Anelastic Liquid and Anelastic Liquid), and a *Timestep Adaptor*, for controlling the time-step length (via a CFL criterion) as the simulation advances in time. "
]
},
{
@@ -304,23 +287,19 @@
"metadata": {},
"outputs": [],
"source": [
- "stress = 2 * mu * sym(grad(u))\n",
- "F_stokes = inner(grad(v), stress) * dx - div(v) * p * dx - (dot(v, k) * Ra * Ttheta) * dx\n",
- "F_stokes += -w * div(u) * dx # Continuity equation\n",
- "F_energy = q * (Tnew - Told) / delta_t * dx + q * dot(u, grad(Ttheta)) * dx + dot(grad(q), kappa * grad(Ttheta)) * dx"
+ "Ra, delta_t = Constant(1e4), Constant(1e-3)\n",
+ "\n",
+ "approximation = BoussinesqApproximation(Ra)\n",
+ "t_adapt = TimestepAdaptor(delta_t, V, maximum_timestep=0.1, increase_tolerance=1.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We specify strong Dirichlet boundary conditions for velocity (bcvx, bcvy) and temperature\n",
- "(bctb, bctt). A Dirichlet boundary condition is created by constructing a Python\n",
- "`DirichletBC` object, where the user must provide the\n",
- "function space the condition applies to, the value, and the part of the mesh at\n",
- "which it applies. \n",
- "Note how boundary conditions are being applied to the velocity part of the\n",
- "mixed finite element space $Z$, indicated by `Z.sub(0)`. Within `Z.sub(0)` we can further subdivide into `Z.sub(0).sub(0)` and `Z.sub(0).sub(1)` to apply boundary conditions to the $x$ and $y$ components of the velocity field only. To apply conditions to the pressure space, we would use `Z.sub(1)`."
+ "We specify strong Dirichlet boundary conditions for velocity and temperature. The user must provide the part of the mesh at\n",
+ "which each boundary condition applies.\n",
+ "Note how boundary conditions have the granularity to be applied to the $x$ and $y$ components of the velocity field only, if desired."
]
},
{
@@ -329,9 +308,17 @@
"metadata": {},
"outputs": [],
"source": [
- "# Set up boundary conditions and deal with nullspaces:\n",
- "bcvx, bcvy = DirichletBC(Z.sub(0).sub(0), 0, sub_domain=(left, right)), DirichletBC(Z.sub(0).sub(1), 0, sub_domain=(bottom, top))\n",
- "bctb, bctt = DirichletBC(Q, 1.0, sub_domain=bottom), DirichletBC(Q, 0.0, sub_domain=top)"
+ "temp_bcs = {\n",
+ " bottom: {\"T\": 1.0},\n",
+ " top: {\"T\": 0.0},\n",
+ "}\n",
+ "\n",
+ "stokes_bcs = {\n",
+ " bottom: {\"uy\": 0},\n",
+ " top: {\"uy\": 0},\n",
+ " left: {\"ux\": 0},\n",
+ " right: {\"ux\": 0},\n",
+ "}"
]
},
{
@@ -340,7 +327,8 @@
"source": [
"With closed boundaries, and no constraint on pressure anywhere in the domain, this problem has a constant pressure nullspace and we must\n",
"ensure that our solver removes this space. To do so, we build a nullspace\n",
- "object, which will subsequently be passed to the solver, and PETSc will seek a solution in the space orthogonal to the provided nullspace. "
+ "object, which will subsequently be passed to the solver, and PETSc will seek a solution in the space orthogonal to the provided nullspace.\n",
+ "When building the nullspace object, the 'closed' keyword handles the constant pressure nullspace, whilst the 'rotational' keyword deals with rotational modes, which, for example, manifest in an a annulus domain with free slip top and bottom boundaries (as we will see in the next tutorial)."
]
},
{
@@ -349,15 +337,15 @@
"metadata": {},
"outputs": [],
"source": [
- "p_nullspace = MixedVectorSpaceBasis(Z, [Z.sub(0), VectorSpaceBasis(constant=True)])"
+ "Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We finally come to solving the variational problem, with associated solver objects. We pass in the residual functions\n",
- "$F_{\\text{Stokes}}$ and $F_{\\text{Energy}}$, solution fields (z, $T_{\\text{new}}$), boundary conditions and, for the Stokes system, the nullspace object. Solution of the two problems is undertaken by PETSc guided by the following solver parameters."
+ "We finally come to solving the variational problem, with solver objects for the energy and Stokes systems created. For the energy system we pass in the solution field T, velocity u, the physical approximation, time step, temporal discretisation approach \n",
+ "(i.e. implicit middle point, being equivalent to the Crank Nicolson scheme outlined above) and boundary conditions. For the Stokes system, we pass in the solution fields z, Temperature, the physical approximation, boundary conditions and the nullspace object. Solution of the two variational problems is undertaken by PETSc. "
]
},
{
@@ -366,34 +354,16 @@
"metadata": {},
"outputs": [],
"source": [
- "# Solver dictionary:\n",
- "solver_parameters = {\n",
- " \"mat_type\": \"aij\",\n",
- " \"snes_type\": \"ksponly\",\n",
- " \"ksp_type\": \"preonly\",\n",
- " \"pc_type\": \"lu\",\n",
- " \"pc_factor_mat_solver_type\": \"mumps\",\n",
- "}\n",
- "\n",
- "# Setup problem and solver objects so we can reuse (cache) solver setup\n",
- "stokes_problem = NonlinearVariationalProblem(F_stokes, z, bcs=[bcvx, bcvy])\n",
- "stokes_solver = NonlinearVariationalSolver(stokes_problem, solver_parameters=solver_parameters, nullspace=p_nullspace, transpose_nullspace=p_nullspace)\n",
- "energy_problem = NonlinearVariationalProblem(F_energy, Tnew, bcs=[bctb, bctt])\n",
- "energy_solver = NonlinearVariationalSolver(energy_problem, solver_parameters=solver_parameters)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "The first option, instructs the Jacobian to be assembled in PETSc's default `aij` sparse matrix type. Although the Stokes and energy problem in this example are linear, for consistency with more realistic cases, we use Firedrake's `NonlinearVariationalSolver` which makes use of PETSc's Scalable Nonlinear Equations Solvers (SNES) interface. However, since we do not actually need a nonlinear solver for this case, we choose the `ksponly` method, indicating that only a single linear solve needs to be performed. The linear solvers are configured through PETSc's Krylov Subspace (KSP) interface, where we can request a direct solver by choosing the `preonly` KSP method, in combination with `lu` as the preconditioner (PC) type. The specific implementation of the LU-decomposition based direct solver is selected as the MUMPS library. Note that the solution process is fully programmable, enabling the creation of sophisticated solvers by combining multiple layers of Krylov methods and preconditioners."
+ "energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)\n",
+ "stokes_solver = StokesSolver(z, T, approximation, bcs=stokes_bcs,\n",
+ " nullspace=Z_nullspace, transpose_nullspace=Z_nullspace)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We can now initiate the time-loop, with the Stokes system solved seperately. These `solve` calls once again convert symbolic mathematics into computation. In the time loop, set here to run for 500 constant time-steps, we compute the RMS velocity and surface Nusselt number for diagnostic purposes, and print these results every 50 timesteps. "
+ "We can now initiate the time-loop, with the Stokes and energy systems solved seperately. These `solve` calls once again convert symbolic mathematics into computation. In the time loop, set here to run for 500 time-steps, we compute the RMS velocity and surface Nusselt number for diagnostic purposes, and print these results every 50 timesteps."
]
},
{
@@ -403,15 +373,18 @@
"outputs": [],
"source": [
"no_timesteps = 500\n",
- "for timestep in range(0, no_timesteps):\n",
+ "time = 0\n",
+ "gd = GeodynamicalDiagnostics(u, p, T, bottom, top)\n",
+ "\n",
+ "for timestep in range(0, no_timesteps+1):\n",
+ " dt = t_adapt.update_timestep(u)\n",
+ " time += dt\n",
+ "\n",
" stokes_solver.solve()\n",
" energy_solver.solve()\n",
- " vrms = sqrt(assemble(dot(u, u) * dx)) * sqrt(1./assemble(1.*dx(domain=mesh)))\n",
- " nu_top = -1. * assemble(dot(grad(Tnew), FacetNormal(mesh)) * ds(top))\n",
+ "\n",
" if timestep % 50 == 0:\n",
- " print(timestep, vrms, nu_top)\n",
- " Told.assign(Tnew)\n",
- "print(timestep, vrms, nu_top)"
+ " print(timestep, dt, gd.u_rms(), gd.Nu_top())"
]
},
{
@@ -425,39 +398,29 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Exercises\n",
- "\n",
- "### 1a: Plot the steady-state temperature field (Tnew), after completion of the 500 timesteps. \n",
- "\n",
- "Do you see a single convective cell? You should!\n",
- "\n",
- "### 1b: Change your initial condition (remembering to keep it bounded between 0 and 1).\n",
- "\n",
- "Do results converge towards the same solution? They should - this is a characteristic of low Rayleigh number Stokes flow. Verify that answers are consistent, even if the upwelling moves to the opposite side of the domain.\n",
- "\n",
- "### 1c: Open this google sheet -- https://docs.google.com/spreadsheets/d/1P45xE0AOv6nkj2HmoXqOJvc0Yrgg4wjU9vSYW7wBLUA/edit?usp=sharing -- pick a resolution, place your name next to it, and run the case at that resolution. \n",
+ "## Exercises\n",
"\n",
- "Add your steady-state Nusselt Number to the relevant column of the spreadsheet. As a group, you should find that results converge towards the benchmark solutions of Blankenbach et al. (1989), with increasing resolution. \n",
+ "1. Plot the steady-state temperature field (T), after completion of the 500 timesteps. Do you see a single convective cell? You should!\n",
"\n",
- "### 1d: Use a lower approximation degree for the temperature field (1) and re-run your cases. What do you notice?\n",
+ "2. Change your initial condition (remembering to keep it bounded between 0 and 1). Do results converge towards the same solution? They should - this is a characteristic of low Rayleigh number Stokes flow. Verify that answers are consistent, even if the upwelling moves to the opposite side of the domain.\n",
"\n",
- "Solve the same problem, only this time, use a piecewise linear approximation space for temperature. How do your results change? Is this as would be expected? \n",
+ "3. Open this google sheet -- https://docs.google.com/spreadsheets/d/1P45xE0AOv6nkj2HmoXqOJvc0Yrgg4wjU9vSYW7wBLUA/edit?usp=sharing -- pick a resolution, place your name next to it, and run the case at that resolution. Add your steady-state Nusselt Number to the relevant column of the spreadsheet. As a group, you should find that results converge towards the benchmark solutions of Blankenbach et al. (1989), with increasing resolution.\n",
"\n",
- "- Hint: check the help for `FunctionSpace` to see how to specify the degree."
+ "4. Solve the same problem, only this time, use a piecewise linear approximation space for temperature (Q1 - hint: check the help for `FunctionSpace` to see how to specify the degree). How do your results change? Is this as would be expected?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Part II: variable viscosity flows"
+ "## Part II: variable viscosity flows"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "In the previous example, we considered isoviscous flows. This is not the case in the mantle, where the viscosity strongly depends on temperature, pressure, composition, strain-rate (non-linear) and potentially many other factors (e.g. grain size). In the second part of this notebook, we look at how straightforward it is to run a case where viscosity depends on temperature, depth and strain-rate -- a so-called visco-plastic rheology -- using a case from Tosi et al. (2015), a benchmark study intended to form a straightforward extension to Blankenbach et al. (1989). Indeed, aside from the viscosity and reference Rayleigh Number ($Ra_{0}=10^2$), all other aspects of this case are identical to the case presented above. \n",
+ "In the previous example, we considered isoviscous flows. This is not the case in the mantle, where the viscosity strongly depends on temperature, pressure, composition, strain-rate (nonlinear) and potentially many other factors (e.g. grain size). In the second part of this notebook, we look at how straightforward it is to run a case where viscosity depends on temperature, depth and strain-rate -- a so-called visco-plastic rheology -- using a case from Tosi et al. (2015), a benchmark study intended to form a straightforward extension to Blankenbach et al. (1989). Indeed, aside from the viscosity and reference Rayleigh Number ($Ra_{0}=10^2$), all other aspects of this case are identical to the case presented above.\n",
"\n",
"The viscosity field, $\\mu$, is calculated as the harmonic mean between a linear component, $\\mu_{\\text{lin}}$ (dependent on depth and temperature) and a nonlinear, plastic component, $\\mu_{\\text{plast}}$, which is dependent on the second invariant of the strain-rate, as follows:\n",
"$$\n",
@@ -471,85 +434,34 @@
"$$\n",
"\\mu_{\\text{plast}}(\\dot\\epsilon) = \\mu^{\\star} + \\frac{\\sigma_{y}}{\\sqrt{\\dot\\epsilon : \\dot\\epsilon}}\n",
"$$\n",
- "where $\\mu^\\star$ is a constant representing the effective viscosity at high stresses and $\\sigma_{y}$ is the yield stress. The denominator of the second term represents the second invariant of the strain-rate tensor. The viscoplastic flow law leads to linear viscous deformation at low stresses and plastic deformation at stresses that exceed $\\sigma_{y}$, with the decrease in viscosity limited by the choice of $\\mu^\\star$. \n",
+ "where $\\mu^\\star$ is a constant representing the effective viscosity at high stresses and $\\sigma_{y}$ is the yield stress. The denominator of the second term represents the second invariant of the strain-rate tensor. The viscoplastic flow law leads to linear viscous deformation at low stresses and plastic deformation at stresses that exceed $\\sigma_{y}$, with the decrease in viscosity limited by the choice of $\\mu^\\star$.\n",
"\n",
- "Although Tosi et al. (2015) examined a number of cases, we focus on one here (Case 4: $Ra_{0}=10^2$, $\\Delta\\mu_T=10^5$, $\\Delta\\mu_{y}=10$ and $\\mu^{\\star}=10^{-3}$), which allows us to demonstrate how a temperature-, depth- and strain-rate dependent viscosity is incorporated within Firedrake. The changes required to simulate this case, relative to our base case are as follows.\n",
- "\n",
- "1. Linear solver options are no longer applicable, given the dependence of viscosity on the flow field, through the strain-rate. Accordingly, the solver dictionary is updated to account for the nonlinear nature of our Stokes system. We exploit the SNES, using a setup based on Newton's method (`snes_type: newtonls`) with a secant line search over the L2-norm of the function (`snes_linesearch_type: l2`). As we target a steady-state solution, an absolute tolerance is specified for our nonlinear solver (`snes_atol: 1e-10`).\n",
- "2. Solver options differ between the (nonlinear) Stokes and (linear) energy systems. As such, a separate solver dictionary is specified for solution of the energy equation. Consistent with our base case, we use a direct solver for solution of the energy equation, based on the MUMPS library. Both of these changes are shown below. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "stokes_solver_parameters = {\n",
- " \"mat_type\": \"aij\",\n",
- " \"snes_type\": \"newtonls\",\n",
- " \"snes_linesearch_type\": \"l2\",\n",
- " \"snes_max_it\": 100,\n",
- " \"snes_atol\": 1e-10,\n",
- " \"ksp_type\": \"preonly\",\n",
- " \"pc_type\": \"lu\",\n",
- " \"pc_factor_mat_solver_type\": \"mumps\",\n",
- "}\n",
- " \n",
- "energy_solver_parameters = {\n",
- " \"mat_type\": \"aij\",\n",
- " \"snes_type\": \"ksponly\",\n",
- " \"ksp_type\": \"preonly\",\n",
- " \"pc_type\": \"lu\",\n",
- " \"pc_factor_mat_solver_type\": \"mumps\",\n",
- "}"
+ "Although Tosi et al. (2015) examined a number of cases, we focus on one here (Case 4: $Ra_{0}=10^2$, $\\Delta\\mu_T=10^5$, $\\Delta\\mu_{y}=10$ and $\\mu^{\\star}=10^{-3}$), which allows us to demonstrate how a temperature-, depth- and strain-rate dependent viscosity is incorporated within Firedrake."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "3. Viscosity is calculated as a function of temperature, depth ($\\mu_{\\text{lin}}$) and strain-rate ($\\mu_{\\text{plast}}$), using the constants specified. Linear and nonlinear components are subsequently combined via a harmonic mean. "
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "Ra = Constant(100.) # Rayleigh number \n",
+ "Viscosity is calculated as a function of temperature, depth (both combined to yield $\\mu_{\\text{lin}}$) and strain-rate ($\\mu_{\\text{plast}}$), using the constants specified. Linear and nonlinear components are subsequently combined via a harmonic mean.\n",
+ "\n",
+ "```python\n",
+ "Ra = Constant(100.) # Rayleigh number\n",
"gamma_T, gamma_Z = Constant(ln(10**5)), Constant(ln(10))\n",
"mu_star, sigma_y = Constant(0.001), Constant(1.0)\n",
- "epsilon = sym(grad(u)) # Strain-rate \n",
- "epsii = sqrt(inner(epsilon,epsilon) + 1e-15) # 2nd invariant (with tolerance to ensure stability)\n",
- "mu_lin = exp(-gamma_T*Tnew + gamma_Z*(1 - X[1]))\n",
+ "epsilon = sym(grad(u)) # Strain-rate\n",
+ "epsii = sqrt(inner(epsilon,epsilon) + 1e-10) # 2nd invariant (with tolerance to ensure stability)\n",
+ "mu_lin = exp(-gamma_T*T + gamma_Z*(1 - X[1]))\n",
"mu_plast = mu_star + (sigma_y / epsii)\n",
- "mu = (2. * mu_lin * mu_plast) / (mu_lin + mu_plast)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "4. Updated solver dictionaries are incorporated into their respective solvers."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "stokes_solver = NonlinearVariationalSolver(stokes_problem, solver_parameters=stokes_solver_parameters, nullspace=p_nullspace, transpose_nullspace=p_nullspace)\n",
- "energy_solver = NonlinearVariationalSolver(energy_problem, solver_parameters=energy_solver_parameters)"
+ "mu = (2. * mu_lin * mu_plast) / (mu_lin + mu_plast)\n",
+ "```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Putting all of this together, the entire script is as follows. "
+ "Putting all of this together, the entire script is as follows."
]
},
{
@@ -558,125 +470,95 @@
"metadata": {},
"outputs": [],
"source": [
- "mesh = UnitSquareMesh(30, 30, quadrilateral=True)\n",
+ "mesh = UnitSquareMesh(40, 40, quadrilateral=True)\n",
"left, right, bottom, top = 1, 2, 3, 4 # Boundary IDs\n",
- "dx = dx(degree=6)\n",
- "V = VectorFunctionSpace(mesh, family=\"CG\", degree=2) # Velocity function space (vector)\n",
- "W = FunctionSpace(mesh, family=\"CG\", degree=1) # Pressure function space (scalar)\n",
- "Q = FunctionSpace(mesh, family=\"CG\", degree=2) # Temperature function space (scalar)\n",
+ "\n",
+ "V = VectorFunctionSpace(mesh, \"CG\", degree=2) # Velocity function space (vector)\n",
+ "W = FunctionSpace(mesh, \"CG\", degree=1) # Pressure function space (scalar)\n",
+ "Q = FunctionSpace(mesh, \"CG\", degree=2) # Temperature function space (scalar)\n",
"Z = MixedFunctionSpace([V, W]) # Mixed function space\n",
"\n",
- "v, w = TestFunctions(Z)\n",
- "q = TestFunction(Q)\n",
"z = Function(Z) # a field over the mixed function space Z.\n",
"u, p = split(z) # returns symbolic UFL expression for u and p\n",
- "Told, Tnew = Function(Q, name=\"OldTemp\"), Function(Q, name=\"NewTemp\")\n",
- "Ttheta = 0.5 * Tnew + 0.5 * Told # Temporal discretisation through Crank-Nicolson\n",
"\n",
+ "T = Function(Q, name=\"Temperature\")\n",
"X = SpatialCoordinate(mesh)\n",
- "Told.interpolate(1.0 - X[1] + 0.05 * cos(pi * X[0]) * sin(pi * X[1]))\n",
- "Tnew.assign(Told)\n",
+ "T.interpolate(1.0 - X[1] + 0.05 * cos(pi * X[0]) * sin(pi * X[1]))\n",
+ "\n",
+ "Ra, delta_t = Constant(100), Constant(1.5e-4)\n",
"\n",
- "Ra, kappa, delta_t = Constant(100), Constant(1.0), Constant(1.5e-4)\n",
- "k = Constant((0, 1)) # Unit vector (in direction opposite to gravity)\n",
+ "approximation = BoussinesqApproximation(Ra)\n",
+ "t_adapt = TimestepAdaptor(delta_t, V, maximum_timestep=0.1, increase_tolerance=1.5)\n",
"\n",
"gamma_T, gamma_Z = Constant(ln(10**5)), Constant(ln(10))\n",
"mu_star, sigma_y = Constant(0.001), Constant(1.0)\n",
- "epsilon = sym(grad(u)) # Strain-rate \n",
- "epsii = sqrt(inner(epsilon,epsilon) + 1e-15) # 2nd invariant (with tolerance to ensure stability)\n",
- "mu_lin = exp(-gamma_T*Tnew + gamma_Z*(1 - X[1]))\n",
+ "epsilon = sym(grad(u)) # Strain-rate\n",
+ "epsii = sqrt(inner(epsilon, epsilon) + 1e-10) # 2nd invariant (with tolerance to ensure stability)\n",
+ "mu_lin = exp(-gamma_T*T + gamma_Z*(1 - X[1]))\n",
"mu_plast = mu_star + (sigma_y / epsii)\n",
"mu = (2. * mu_lin * mu_plast) / (mu_lin + mu_plast)\n",
"\n",
- "stress = 2 * mu * sym(grad(u))\n",
- "F_stokes = inner(grad(v), stress) * dx - div(v) * p * dx - (dot(v, k) * Ra * Ttheta) * dx\n",
- "F_stokes += -w * div(u) * dx # Continuity equation\n",
- "F_energy = q * (Tnew - Told) / delta_t * dx + q * dot(u, grad(Ttheta)) * dx + dot(grad(q), kappa * grad(Ttheta)) * dx\n",
- "\n",
- "bcvx, bcvy = DirichletBC(Z.sub(0).sub(0), 0, sub_domain=(left, right)), DirichletBC(Z.sub(0).sub(1), 0, sub_domain=(bottom, top))\n",
- "bctb, bctt = DirichletBC(Q, 1.0, sub_domain=bottom), DirichletBC(Q, 0.0, sub_domain=top)\n",
- "p_nullspace = MixedVectorSpaceBasis(Z, [Z.sub(0), VectorSpaceBasis(constant=True)])\n",
- "\n",
- "stokes_solver_parameters = {\n",
- " \"mat_type\": \"aij\",\n",
- " \"snes_type\": \"newtonls\",\n",
- " \"snes_linesearch_type\": \"l2\",\n",
- " \"snes_max_it\": 100,\n",
- " \"snes_atol\": 1e-10,\n",
- " \"ksp_type\": \"preonly\",\n",
- " \"pc_type\": \"lu\",\n",
- " \"pc_factor_mat_solver_type\": \"mumps\",\n",
+ "Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=False)\n",
+ "\n",
+ "temp_bcs = {\n",
+ " bottom: {\"T\": 1.0},\n",
+ " top: {\"T\": 0.0},\n",
"}\n",
- " \n",
- "energy_solver_parameters = {\n",
- " \"mat_type\": \"aij\",\n",
- " \"snes_type\": \"ksponly\",\n",
- " \"ksp_type\": \"preonly\",\n",
- " \"pc_type\": \"lu\",\n",
- " \"pc_factor_mat_solver_type\": \"mumps\",\n",
+ "\n",
+ "stokes_bcs = {\n",
+ " bottom: {\"uy\": 0},\n",
+ " top: {\"uy\": 0},\n",
+ " left: {\"ux\": 0},\n",
+ " right: {\"ux\": 0},\n",
"}\n",
"\n",
- "stokes_problem = NonlinearVariationalProblem(F_stokes, z, bcs=[bcvx, bcvy])\n",
- "stokes_solver = NonlinearVariationalSolver(stokes_problem, solver_parameters=stokes_solver_parameters, nullspace=p_nullspace, transpose_nullspace=p_nullspace)\n",
- "energy_problem = NonlinearVariationalProblem(F_energy, Tnew, bcs=[bctb, bctt])\n",
- "energy_solver = NonlinearVariationalSolver(energy_problem, solver_parameters=energy_solver_parameters)\n",
+ "energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)\n",
+ "stokes_solver = StokesSolver(z, T, approximation, bcs=stokes_bcs, mu=mu,\n",
+ " nullspace=Z_nullspace, transpose_nullspace=Z_nullspace)\n",
"\n",
"no_timesteps = 2000\n",
- "for timestep in range(0, no_timesteps):\n",
+ "time = 0\n",
+ "gd = GeodynamicalDiagnostics(u, p, T, bottom, top)\n",
+ "\n",
+ "for timestep in range(0, no_timesteps+1):\n",
+ " dt = t_adapt.update_timestep(u)\n",
+ " time += dt\n",
+ "\n",
" stokes_solver.solve()\n",
" energy_solver.solve()\n",
- " vrms = sqrt(assemble(dot(u, u) * dx)) * sqrt(1./assemble(1.*dx(domain=mesh)))\n",
- " nu_top = -1. * assemble(dot(grad(Tnew), FacetNormal(mesh)) * ds(top))\n",
+ "\n",
" if timestep % 50 == 0:\n",
- " print(timestep, vrms, nu_top)\n",
- " Told.assign(Tnew)\n",
- "print(timestep, vrms, nu_top)"
+ " print(timestep, dt, gd.u_rms(), gd.Nu_top())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Exercises\n",
- "\n",
- "### 2a: Plot the steady-state viscosity field. \n",
- "\n",
- "Note that you will need to create a field and interpolate values to that field, as the viscosity $\\mu$ is currently only specified in UFL.\n",
+ "## Exercises\n",
"\n",
- "### 2b: Open this google sheet -- https://docs.google.com/spreadsheets/d/1P45xE0AOv6nkj2HmoXqOJvc0Yrgg4wjU9vSYW7wBLUA/edit?usp=sharing -- navigate to the \"Tosi_Case\" tab, pick a resolution, place your name next to it, and run the case at that resolution. \n",
+ "5. Plot the steady-state viscosity field. Note that you will need to create a field and interpolate values on to that field, as the viscosity $\\mu$ is currently only specified in UFL.\n",
"\n",
- "Add your steady-state Nusselt Number to the relevant column of the spreadsheet. As a group, you should find that results converge towards the benchmark solutions of Tosi et al. (2015), with increasing resolution. \n",
+ "6. Open this google sheet -- https://docs.google.com/spreadsheets/d/1P45xE0AOv6nkj2HmoXqOJvc0Yrgg4wjU9vSYW7wBLUA/edit?usp=sharing -- navigate to the \"Tosi_Case\" tab, pick a resolution, place your name next to it, and run the case at that resolution. Add your steady-state Nusselt Number to the relevant column. As a group, you should find that results converge towards the benchmark solutions of Tosi et al. (2015), with increasing resolution.\n",
"\n",
- "### 2c: Increase the sensitivity of viscosity to temperature. \n",
- "\n",
- "What do you observe? \n",
- "\n"
+ "7. Increase the sensitivity of viscosity to temperature. What do you observe?"
]
}
],
"metadata": {
"interpreter": {
- "hash": "2834c7709b7175be062475c9c8d239d84f48b92a4db284519cf4bc00989e411b"
+ "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
},
"kernelspec": {
- "display_name": "Python 3.9.12 ('firedrake')",
+ "display_name": "Python 3.11.5 64-bit",
"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.13"
- },
- "orig_nbformat": 4
+ "version": "3.11.5"
+ }
},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/07-GD-2D-cylindrical.ipynb b/07-GD-2D-cylindrical.ipynb
index 32092e4..9eb13af 100644
--- a/07-GD-2D-cylindrical.ipynb
+++ b/07-GD-2D-cylindrical.ipynb
@@ -1,14 +1,5 @@
{
"cells": [
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "!sudo mv /etc/apt/sources.list.d/cuda.list{,.disabled}"
- ]
- },
{
"cell_type": "code",
"execution_count": null,
@@ -16,26 +7,27 @@
"outputs": [],
"source": [
"# This magic makes plots appear \n",
- "%matplotlib notebook\n",
"%matplotlib inline\n",
- "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "import os\n",
+ "os.environ[\"GADOPT_LOGLEVEL\"] = \"WARNING\"\n",
+ "\n",
"# Load Firedrake on Colab\n",
"try:\n",
" import firedrake\n",
"except ImportError:\n",
- " !wget \"https://fem-on-colab.github.io/releases/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
- " import firedrake\n",
- "from firedrake.petsc import PETSc\n",
- "import matplotlib.pyplot as plt"
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " import firedrake"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "# 2-D Cylindrical Case\n",
- "## Weak formulation \n",
- "For our reference we take another look at our weak formulation: \n",
+ "# 2-D Annulus Case\n",
+ "## Weak formulation\n",
+ "As outlined in the previous tutorial, the weak formulation of the equations governing mantle dynamics are: \n",
"\n",
"$$\\int_\\Omega (\\nabla \\vec{v})\\colon \\mu \\left[ \\nabla \\vec{u} + \\left( \\nabla \\vec{u} \\right)^T\\right] \\ dx \n",
" - \\int_{\\Omega} \\left( \\nabla \\cdot \\vec{v}\\right)\\ p \\ dx\n",
@@ -47,33 +39,33 @@
" + \\int_\\Omega q \\vec{u}\\cdot \\nabla T \\ dx \n",
" + \\int_\\Omega \\left(\\nabla q\\right) \\cdot \\left(\\kappa \\nabla T\\right) \\ dx = 0 \\text{ for all } q\\in Q.$$\n",
"\n",
+ "This example analyses mantle flow, subject to these same governing equations, in a 2-D annulus domain. We define our domain by the radii of the inner ($r_{\\text{min}}$) and outer ($r_{\\text{max}}$) boundaries. \n",
+ "These are chosen such that the non--dimensional depth of the mantle, $z = r_{\\text{max}} - r_{\\text{min}} = 1$, and the ratio of the inner and outer radii, $f=r_{\\text{min}} / r_{\\text{max}} = 0.55$, \n",
+ "thus approximating the ratio between the radii of Earth's surface and core-mantle-boundary (CMB). Specifically, we set $r_{\\text{min}} = 1.22$ and $r_{\\text{max}} = 2.22$. \n",
"\n",
- "In the Cartesian examples we saw earlier are imposed through strong Dirichlet boundary conditions for velocity $\\vec{u}$. This is achieved by restricting the velocity function space $V$ to a subspace $V_0$\n",
- "of vector functions for which all components (zero-slip) or only the normal\n",
- "component (free-slip) are zero at the boundary. Since this restriction also\n",
- "applies to the test functions $\\vec{v}$, the weak form only needs to be\n",
- "satisfied for all test functions $\\vec{v}\\in V_0$ that satisfy the homogeneous boundary conditions. Therefore, the omitted boundary integral\n",
- "\\begin{equation}\n",
- " -\\int_{\\partial\\Omega} \\vec{v}\\cdot \\left(\\mu \\left[\\nabla\\vec{u}\n",
- " + \\left(\\nabla\\vec{u}\\right)^T\\right]\\right)\\cdot \\vec{n} ds\n",
- "\\end{equation}\n",
- "that was required to obtain the integrated by parts viscosity term, automatically vanishes for zero-slip boundary conditions as\n",
- "$\\bf v =0$ at the domain boundary, $\\partial\\Omega$. In the case of a free-slip\n",
- "boundary condition for which the tangential components of $\\vec{v}$ are\n",
- "non-zero, the boundary term does not vanish, but by omitting that term, we weakly impose a zero shear stress condition. The boundary\n",
- "term obtained by integrating the pressure gradient term by parts,\n",
- "\\begin{equation}\n",
- " \\int_{\\partial\\Omega} \\vec{v}\\cdot\\vec{n} p ds ,\n",
- "\\end{equation}\n",
- "\n",
- "\n",
- "Note that we have integrated by parts the viscosity and pressure gradient terms in the Stokes equations, and the diffusion term in the energy equation, but have omitted the corresponding boundary terms. Our goal in this section is to examine simulations in a 2-D cylindrical domain. \n",
- "\n",
- "We define our domain by the radii of the inner ($r_{\\text{min}}$) and outer ($r_{\\text{max}}$) boundaries. These are chosen such that the non--dimensional depth of the mantle,\n",
- "$z = r_{\\text{max}} - r_{\\text{min}} = 1$, and the ratio of the inner and outer radii, $f=r_{\\text{min}} / r_{\\text{max}} = 0.55$, thus approximating the ratio between the\n",
- "radii of Earth's surface and core-mantle-boundary (CMB). Specifically, we set $r_{\\text{min}} = 1.22$ and $r_{\\text{max}} = 2.22$. \n",
- "\n",
- "\n"
+ "This example focusses on differences between running simulations in a 2-D annulus and 2-D Cartesian domain. These can be summarised as follows:\n",
+ "1. The geometry of the problem - i.e. the computational mesh. \n",
+ "2. The radial direction of gravity (as opposed to the vertical direction in a Cartesian domain).\n",
+ "3. Initialisation of the temperature field in a domain of different dimensions.\n",
+ "4. With free-slip boundary conditions on both boundaries, this case incorporates a velocity nullspace, as well as a pressure nullspace.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from gadopt import *\n",
+ "import numpy as np\n",
+ "import pyvista as pv # Used for plotting vtk output"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We first need a mesh. We generate a circular manifold mesh (with 32 elements) and extrude in the radial direction, using the optional keyword argument `extrusion_type`, forming 8 layers. To better represent the curvature of the domain and ensure accuracy of our quadratic representation of velocity, we approximate the curved cylindrical shell domain quadratically, using the optional keyword argument `degree`$=2$."
]
},
{
@@ -82,49 +74,22 @@
"metadata": {},
"outputs": [],
"source": [
- "from firedrake import *\n",
- "from mpi4py import MPI\n",
- "# Quadrature degree for our integrals:\n",
- "dx = dx(degree=6)\n",
- "\n",
"# Set up geometry:\n",
"rmin, rmax, ncells, nlayers = 1.22, 2.22, 32, 8\n",
"\n",
- "# Function to generate a cylindrical mesh\n",
- "# cylinder_mesh is replacing the following \n",
- "#mesh1d = CircleManifoldMesh(ncells, radius=rmin, degree=2)\n",
- "#mesh = ExtrudedMesh(mesh1d, layers=nlayers, extrusion_type='radial')\n",
- "def cylinder_mesh(rmin, rmax, ncells, nlayers):\n",
- " cmesh = CylinderMesh(ncells, nlayers, quadrilateral=True)\n",
- " coord_fs = VectorFunctionSpace(cmesh, FiniteElement(\"Q\", 'quadrilateral', 1, variant=\"equispaced\"), dim=2)\n",
- " new_coords = Function(coord_fs)\n",
- " new_coords.dat.data[:] = cmesh.coordinates.dat.data[:,:2] * (rmin + cmesh.coordinates.dat.data[:,2][:,np.newaxis] * (rmax-rmin))\n",
- " mesh = Mesh(new_coords)\n",
- " new_coords2 = Function(VectorFunctionSpace(mesh, \"CG\", 2))\n",
- " new_coords2.interpolate(SpatialCoordinate(mesh))\n",
- " x, y = SpatialCoordinate(mesh)\n",
- " rlin = interpolate(sqrt(x**2 + y**2), FunctionSpace(mesh, \"CG\", 1))\n",
- " new_coords2.interpolate(new_coords2/sqrt(new_coords2[0]**2 + new_coords2[1]**2) * rlin)\n",
- " return Mesh(new_coords2)\n",
+ "mesh1d = CircleManifoldMesh(ncells, radius=rmin, degree=2)\n",
+ "mesh = ExtrudedMesh(mesh1d, layers=nlayers, extrusion_type='radial')\n",
"\n",
- "\n",
- "# generating mesh\n",
- "mesh = cylinder_mesh(rmin, rmax, ncells, nlayers)\n",
- "bottom_id, top_id = 1, 2 # the ids for the bottom surface and the top one \n",
- "n = FacetNormal(mesh) # Normals, required for Nusselt number calculation\n",
- "domain_volume = assemble(1*dx(domain=mesh)) # Required for diagnostics (e.g. RMS velocity)\n",
- "top_surface = assemble(Constant(1, domain=mesh)*ds(1))\n",
- "bottom_surface = assemble(Constant(1, domain=mesh)*ds(2))"
+ "bottom_id, top_id = \"bottom\", \"top\"\n",
+ "n = FacetNormal(mesh) # used later for diagnostic purposes\n",
+ "domain_volume = assemble(1 * dx(domain=mesh)) # used later for diagnostic purposes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Exercise 7.1:\n",
- "- Visualise the generated mesh.\n",
- "\n",
- "- Change radial and lateral resolution and visualise the mesh again. \n"
+ "We can now visualise the resulting mesh:"
]
},
{
@@ -133,16 +98,32 @@
"metadata": {},
"outputs": [],
"source": [
- "fig, axes = plt.subplots()\n",
- "triplot(mesh, axes=axes)\n",
- "axes.legend();"
+ "V = FunctionSpace(mesh, \"CG\", 1)\n",
+ "File(\"mesh.pvd\").write(Function(V))\n",
+ "mesh_data = pv.read(\"mesh_0.vtu\")\n",
+ "edges = mesh_data.extract_all_edges()\n",
+ "\n",
+ "plotter = pv.Plotter(notebook=True)\n",
+ "plotter.add_mesh(edges, color=\"black\")\n",
+ "plotter.camera_position = \"xy\"\n",
+ "\n",
+ "plotter.show(jupyter_backend=\"static\", interactive=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "We now need the function spaces that are associated with this mesh. We follow the same definitions as our last exercise, and use the following definitions: \n"
+ "### Exercise:\n",
+ "\n",
+ "1. Change the number of radial layers (nlayers) and the lateral resolution (ncells) and visualise the mesh again."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now need the function spaces associated with this mesh, test functions and functions to hold our solutions. We follow the same definitions as our previous tutorial, using the bilinear Q2-Q1 element pair for velocity and pressure and a Q2 discretisation for temperature."
]
},
{
@@ -151,16 +132,15 @@
"metadata": {},
"outputs": [],
"source": [
- "# Set up function spaces - currently using the bilinear Q2Q1 element pair:\n",
- "V = VectorFunctionSpace(mesh, \"CG\", 2) # Velocity function space (vector)\n",
- "W = FunctionSpace(mesh, \"CG\", 1) # Pressure function space (scalar)\n",
- "Q = FunctionSpace(mesh, \"CG\", 2) # Temperature function space (scalar)\n",
- "Z = MixedFunctionSpace([V, W]) # Mixed function space.\n",
- "# Test functions and functions to hold solutions:\n",
+ "V = VectorFunctionSpace(mesh, \"CG\", degree=2) # Velocity function space (vector)\n",
+ "W = FunctionSpace(mesh, \"CG\", degree=1) # Pressure function space (scalar)\n",
+ "Q = FunctionSpace(mesh, \"CG\", degree=2) # Temperature function space (scalar)\n",
+ "Z = MixedFunctionSpace([V, W]) # Mixed function space\n",
"v, w = TestFunctions(Z)\n",
"q = TestFunction(Q)\n",
"z = Function(Z) # a field over the mixed function space Z.\n",
- "u, p = split(z) # Returns symbolic UFL expression for u and p"
+ "u, p = split(z) # returns symbolic UFL expression for u and p\n",
+ "T = Function(Q, name=\"Temperature\")"
]
},
{
@@ -168,8 +148,9 @@
"metadata": {},
"source": [
"\n",
- "We choose the initial temperature distribution in a way so that our convection run produces 4 equidistant plumes. This initial temperature field is prescribed as:\n",
+ "We choose the initial temperature distribution to trigger upwelling of 4 equidistant plumes. This initial temperature field is prescribed as:\n",
"$$T(x,y) = (r_{\\text{max}} - r) + A\\cos(4 \\; atan2\\ (y,x)) \\sin(r-r_{\\text{min}}) \\pi)$$\n",
+ "\n",
"where $A=0.02$ is the amplitude of the initial perturbation. "
]
},
@@ -182,30 +163,42 @@
"# Set up temperature field and initialise:\n",
"X = SpatialCoordinate(mesh)\n",
"r = sqrt(X[0]**2 + X[1]**2)\n",
- "Told, Tnew = Function(Q, name=\"OldTemp\"), Function(Q, name=\"NewTemp\");\n",
- "pi = 3.141592653589793238\n",
- "Told.interpolate(rmax - r + 0.02*cos(4*atan_2(X[1], X[0])) * sin((r - rmin) * pi))\n",
- "Tnew.assign(Told)"
+ "T.interpolate(rmax - r + 0.02*cos(4*atan_2(X[1], X[0])) * sin((r - rmin) * pi))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can now visualise this initial temperature field:"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# Let's visualise the initial temperature field now\n",
- "# NBVAL_IGNORE_OUTPUT\n",
- "fig, axes = plt.subplots()\n",
- "collection = tripcolor(Tnew, axes=axes, cmap='coolwarm')\n",
- "fig.colorbar(collection);"
+ "File(\"temp.pvd\").write(T)\n",
+ "temp_data = pv.read(\"temp_0.vtu\")\n",
+ "plotter = pv.Plotter(notebook=True)\n",
+ "plotter.add_mesh(temp_data)\n",
+ "plotter.camera_position = \"xy\"\n",
+ "\n",
+ "plotter.show(jupyter_backend=\"static\", interactive=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Next, we strategise our solutions for the energy equations. We use $\\theta$ time-stepping scheme for our enery solutions, and define our way of computing a time step consistent with a cfl condition number of 1.0."
+ "The Rayleigh number for this problem is defined in addition to the initial timestep $\\Delta t$. The viscosity and thermal diffusivity are left at their default values (both = 1). \n",
+ "These constants are used to create an *Approximation* representing the physical setup of the problem (options include Boussinesq, Extended Boussinesq, Truncated Anelastic Liquid \n",
+ "and Anelastic Liquid), and a *Timestep Adaptor*, for controlling the time-step length (via a CFL criterion) as the simulation advances in time. "
]
},
{
@@ -214,46 +207,19 @@
"metadata": {},
"outputs": [],
"source": [
- "# Temporal discretisation - Using a Crank-Nicholson scheme where theta = 0.5:\n",
- "Ttheta = 0.5*Tnew + (1 - 0.5)*Told\n",
- "\n",
- "# Define time stepping parameters:\n",
- "steady_state_tolerance = 1e-9\n",
- "max_timesteps = 2000 \n",
- "target_cfl_no = 1.0\n",
- "maximum_timestep = 0.1\n",
- "increase_tolerance = 1.5\n",
- "time = 0.0\n",
- "\n",
- "# Timestepping - CFL related stuff:\n",
- "ref_vel = Function(V, name=\"Reference_Velocity\")\n",
+ "Ra, delta_t = Constant(1e5), Constant(1e-6)\n",
"\n",
- "def compute_timestep(u, current_delta_t):\n",
- " \"\"\"Return the timestep, based upon the CFL criterion\"\"\"\n",
- "\n",
- " ref_vel.interpolate(dot(JacobianInverse(mesh), u))\n",
- " ts_min = 1. / mesh.comm.allreduce(ref_vel.dat.data.max(), MPI.MAX)\n",
- " # Grab (smallest) maximum permitted on all cores:\n",
- " ts_max = min(float(current_delta_t) * increase_tolerance, maximum_timestep)\n",
- " # Compute timestep:\n",
- " tstep = min(ts_min * target_cfl_no, ts_max)\n",
- " return tstep"
+ "approximation = BoussinesqApproximation(Ra)\n",
+ "t_adapt = TimestepAdaptor(delta_t, V, maximum_timestep=0.1, increase_tolerance=1.5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Setting Solve Parameters - An iterative approach\n",
- "Although direct solves provide robust performance in small problems, in 3-D and large problems, the computational (CPU and memory) requirements quickly become intractable. PETSc's `fieldsplit pc_type` provides a class of preconditioners for mixed problems that allows one to apply different preconditioners to different blocks of the system. Here we configure the Schur complement approach (see GMD paper for more info) \n",
- " \n",
- "The `fieldsplit_0` entries configure solver options for the first of these blocks, the $K$ matrix. The linear systems associated with this matrix are solved using a combination of the Conjugate Gradient method (`cg`, line 10) and an algebraic multigrid preconditioner (`gamg`, line 15). We also specify two options (`gamg_threshold` and `gamg_square_graph`) that control the aggregation method (coarsening strategy) in the GAMG preconditioner, which balance the multigrid effectiveness (convergence rate) with coarse grid complexity (cost per iteration). \n",
- "\n",
- "`fieldsplit_1` contains solver options for the Schur complement solve. Here we approxmiate the Schur complement matrix with a mass matrix scaled by viscosity, which is implementd in `MassInvPC`, with the viscosity provided through the option `apptcx` (Line 29 of a few cells below).\n",
- "\n",
- " Specification of the matrix type `matfree` (line 3) for the combined system ensures that we do not explicitly assemble its associated sparse matrix, instead computing the matrix-vector multiplications required by the Krylov iterations as they arise. Again, for preconditioning in the $K$-matrix solve we need access to matrix values, which is achieved using `AssembledPC`. This explicitly assembles the $K$-matrix by extracting relevant terms from the `F_Stokes` form.\n",
- "\n",
- "Finally, the energy solve is performed through a combination of the GMRES (`gmres`) Krylov method and SSOR preconditioning (lines 27). For all iterative solves we specify a convergence criterion based on the relative reduction of the preconditioned residual (`ksp_rtol`)."
+ "Boundary conditions for temperature are set to $T = 0$ at the surface ($r_{\\text{max}}$) and $T = 1$ at the base ($r_{\\text{min}}$).\n",
+ "For velocity, we specify free‐slip conditions on both boundaries. We incorporate these weakly through the Nitsche approximation.\n",
+ "This illustrates a key advantage of the G-ADOPT framework: the user only specifies that the normal component of velocity is zero and all required changes are handled under the hood. "
]
},
{
@@ -262,44 +228,14 @@
"metadata": {},
"outputs": [],
"source": [
- "# Stokes Equation Solver Parameters:\n",
- "stokes_solver_parameters = {\n",
- " \"mat_type\": \"matfree\",\n",
- " \"snes_type\": \"ksponly\",\n",
- " \"ksp_type\": \"preonly\",\n",
- " \"pc_type\": \"fieldsplit\",\n",
- " \"pc_fieldsplit_type\": \"schur\",\n",
- " \"pc_fieldsplit_schur_type\": \"full\",\n",
- " \"fieldsplit_0\": {\n",
- " \"ksp_type\": \"cg\",\n",
- " \"ksp_rtol\": 1e-5,\n",
- " #\"ksp_converged_reason\": None,\n",
- " \"pc_type\": \"python\",\n",
- " \"pc_python_type\": \"firedrake.AssembledPC\",\n",
- " \"assembled_pc_type\": \"gamg\",\n",
- " \"assembled_pc_gamg_threshold\": 0.01,\n",
- " \"assembled_pc_gamg_square_graph\": 100,\n",
- " },\n",
- " \"fieldsplit_1\": {\n",
- " \"ksp_type\": \"fgmres\",\n",
- " \"ksp_rtol\": 1e-4,\n",
- " #\"ksp_converged_reason\": None,\n",
- " \"pc_type\": \"python\",\n",
- " \"pc_python_type\": \"firedrake.MassInvPC\",\n",
- " \"Mp_ksp_rtol\": 1e-5,\n",
- " \"Mp_ksp_type\": \"cg\",\n",
- " \"Mp_pc_type\": \"sor\",\n",
- " }\n",
+ "temp_bcs = {\n",
+ " bottom_id: {\"T\": 1.0},\n",
+ " top_id: {\"T\": 0.0},\n",
"}\n",
"\n",
- "# Energy Equation Solver Parameters:\n",
- "energy_solver_parameters = {\n",
- " \"mat_type\": \"aij\",\n",
- " \"snes_type\": \"ksponly\",\n",
- " \"ksp_type\": \"gmres\",\n",
- " \"ksp_rtol\": 1e-5,\n",
- " #\"ksp_converged_reason\": None,\n",
- " \"pc_type\": \"sor\",\n",
+ "stokes_bcs = {\n",
+ " bottom_id: {\"un\": 0},\n",
+ " top_id: {\"un\": 0},\n",
"}"
]
},
@@ -307,45 +243,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Boundary Condition: Introducting Weak Imposition of Boundary Conditions\n",
- "Boundary conditions for temperature are $T = 0$ at the surface ($r_{\\text{max}}$) and $T = 1$ at the base ($r_{\\text{min}}$). \n",
- "For velocity boundary conditions we choose free‐slip velocity boundary conditions that are specified on both boundaries. We incorporate these weakly through the Nitsche approximation.\n",
- "### Weak Imposition of Boudary Conditions - Nitsche Penalty Method\n",
- "In curved domains, such as the 2-D cylindrical and 3-D spherical cases examined below, imposing free-slip boundary conditions \n",
- "is complicated by the fact that it is not straightforward to decompose the degrees of freedom of the velocity space V\n",
- "into tangential and lateral components for many finite element discretisations. For Lagrangian based discretisations we could\n",
- "define normal vectors at the Lagrangian nodes on the surface and decompose accordingly, but these normal vectors would have\n",
- "to be averaged due to the piecewise approximation of the curved surface. To avoid such complications for our examples in\n",
- "cylindrical and spherical geometries, we employ a symmetric Nitsche penalty method where the velocity space\n",
- "is not restricted and, thus, retains all discrete solutions with a non-zero normal component. \n",
- "\n",
- "This entails adding the following three surface integrals:\n",
- "\\begin{equation}\n",
- " - \\int_{\\partial\\Omega} \\vec{v}\\cdot\\vec{n}~\n",
- " \\vec{n}\\cdot \\left(\\mu \\left[\\nabla\\vec{u}\n",
- " + \\left(\\nabla\\vec{u}\\right)^T\\right]\\right)\\cdot \\vec{n} \\ ds\n",
- " - \\int_{\\partial\\Omega} \n",
- " \\vec{n}\\cdot \\left(\\mu \\left[\\nabla\\vec{v}\n",
- " + \\left(\\nabla\\vec{v}\\right)^T\\right]\\right)\\cdot \\vec{n}~\n",
- " \\vec{u}\\cdot\\vec{n} \\ ds\n",
- " + \\int_{\\partial\\Omega} C_{\\text{Nitsche}} \\, \\mu \\vec{v}\\cdot\\vec{n} ~ \\vec\n",
- " u\\cdot\\vec{n} \\ ds .\n",
- "\\end{equation}\n",
- "\n",
+ "As noted above, with a free-slip boundary condition on both boundaries, one can add an arbitrary rotation of the form $(-y, x)=r\\hat{\\mathbf{\\theta}}$ to the velocity solution (i.e. this case incorporates a velocity nullspace, as well as a pressure nullspace). These lead to null-modes (eigenvectors) for the linear system, rendering the resulting matrix singular. In preconditioned Krylov methods these null-modes must be subtracted from the approximate solution at every iteration. We do that below, setting up a nullspace object as we did in the previous tutorial, albeit speciying the `rotational` keyword argument to be True. Once again, this removes the requirement for a user to configure these options, further simplifying the task of setting up a (valid) geodynamical simulation. \n",
"\n",
- "The first of these corresponds to the normal component of stress \n",
- "associated with integration by parts of the viscosity term. The tangential\n",
- "component, as before, is omitted and weakly imposes a zero shear stress\n",
- "condition. The second term ensures symmetry of Equation \\eqref{eq:weak_mom} \n",
- "with respect to $\\vec{u}$ and $\\vec{v}$. The third term penalizes the normal\n",
- "component of $\\vec{u}$ and involves a penalty parameter $C_{\\text{Nitsche}}>0$\n",
- "that should be sufficiently large to ensure convergence.\n",
- "\n",
- " \n",
- "coercivity of the bilinear form $F_{\\text{Stokes}}$ introduced in\n",
- "Lower bounds for $C_{\\text{Nitsche},f}$ on each face $f$ can be derived for simplicial \\citep{shahbazi_2005} and quadrilateral/hexahedral \\citep{hillewaert_2013} meshes\n",
- "\n",
- "Finally for the energy equation, we apply a simple $\\theta$ scheme "
+ "Given the increased computational expense (typically requiring more degrees of freedom) in a 2-D annulus domain, the G-ADOPT library defaults to iterative solver parameters (in 2-D Cartesian domains, the framework defaults to direct solver parameters). Our iterative solver setup is configured to use the GAMG preconditioner for the velocity block of the Stokes system, to which we must provide near-nullspace information, which consists of three rotational (`x_rotV`, `y_rotV`, `z_rotV`) and three translational (`nns_x`, `nns_y`, `nns_z`) modes (see Davies et al. GMD, 2022 for details). "
]
},
{
@@ -354,48 +254,17 @@
"metadata": {},
"outputs": [],
"source": [
- "\n",
- "# Stokes related constants (note that since these are included in UFL, they are wrapped inside Constant):\n",
- "Ra = Constant(1e5) # Rayleigh number\n",
- "mu = Constant(1.0) # Viscosity\n",
- "k = as_vector((X[0], X[1])) / r\n",
- "C_ip = Constant(100.0) # Fudge factor for interior penalty term used in weak imposition of BCs\n",
- "p_ip = 2 # Maximum polynomial degree of the _gradient_ of velocity\n",
- "\n",
- "# Temperature equation related constants:\n",
- "delta_t = Constant(1e-7) # Initial time-step\n",
- "kappa = Constant(1.0) # Thermal diffusivity\n",
- "\n",
- "# Stokes equations in UFL form:\n",
- "def stress(u):\n",
- " return 2 * mu * sym(grad(u))\n",
- "# EXCERSIE 7.2 & 7.3\n",
- "F_stokes = inner(grad(v), stress(u)) * dx - div(v) * p * dx + dot(n, v) * p * (ds(top_id)+ds(bottom_id)) - (dot(v, k) * Ra * Ttheta) * dx\n",
- "F_stokes += -w * div(u) * dx + w * dot(n, u) * (ds(top_id)+ds(bottom_id)) # Continuity equation\n",
- "\n",
- "# nitsche free-slip BCs\n",
- "F_stokes += -dot(v, n) * dot(dot(n, stress(u)), n) * (ds(top_id)+ds(bottom_id))\n",
- "F_stokes += -dot(u, n) * dot(dot(n, stress(v)), n) * (ds(top_id)+ds(bottom_id))\n",
- "F_stokes += C_ip * mu * (p_ip + 1)**2 * FacetArea(mesh) / CellVolume(mesh) * dot(u, n) * dot(v, n) * (ds(top_id)+ds(bottom_id))\n",
- "\n",
- "# Energy equation in UFL form:\n",
- "F_energy = q * (Tnew - Told) / delta_t * dx + q * dot(u, grad(Ttheta)) * dx + dot(grad(q), kappa * grad(Ttheta)) * dx\n",
- "\n",
- "# Temperature boundary conditions\n",
- "bctb, bctt = DirichletBC(Q, 1.0, bottom_id), DirichletBC(Q, 0.0, top_id)\n",
- "#bcvt = DirichletBC(Z.sub(0), Constant((0,0)), top_id)\n"
+ "Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=True)\n",
+ "Z_near_nullspace = create_stokes_nullspace(Z, closed=False, rotational=True, translations=[0, 1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Nullspaces and near-nullspaces \n",
- "With a free-slip boundary condition on both boundaries, one can add an arbitrary rotation of the form $(-y, x)=r\\hat{\\mathbf{\\theta}}$ to the velocity solution (i.e. this case incorporates a velocity nullspace, as well as a pressure nullspace). These lead to null-modes (eigenvectors) for the linear system, rendering the resulting matrix singular. In preconditioned Krylov methods these null-modes must be subtracted from the approximate solution at every iteration.\n",
- "\n",
- "By using the GAMG preconditioner, we must provide near-nullspace information for the GAMG preconditioner, consisting of three rotational (`x_rotV`, `y_rotV`, `z_rotV`) and three translational (`nns_x`, `nns_y`, `nns_z`) modes (see GMD paper for details). \n",
- "\n",
- "We now update the Stokes problem to account for additional boundary conditions, and the Stokes solver to include the near nullspace options defined above, in addition to the optional `appctx` keyword argument that passes the viscosity through to our `MassInvPC` Schur complement preconditioner. Energy solver options are also updated relative to our base case, using the dictionary that is created before.\n"
+ "We next come to solving the variational problem, with solver objects for the energy and Stokes systems created. For the energy system we pass in the solution field T, velocity u, the physical approximation, time step, temporal discretisation approach \n",
+ "(i.e. implicit middle point, being equivalent to a Crank Nicolson scheme) and boundary conditions. For the Stokes system, we pass in the solution fields z, Temperature, the physical approximation, boundary condition, the nullspace and near nullspace objects. We also set the\n",
+ "optional `cartesian` keyword argument to False, which ensures that the unit vector, $\\hat{k}$, points radially, in the direction opposite to gravity. Solution of the two variational problems is undertaken by PETSc. "
]
},
{
@@ -404,48 +273,18 @@
"metadata": {},
"outputs": [],
"source": [
- "\n",
- "# Nullspaces and near-nullspaces:\n",
- "x_rotV = Function(V).interpolate(as_vector((-X[1], X[0])))\n",
- "V_nullspace = VectorSpaceBasis([x_rotV])\n",
- "V_nullspace.orthonormalize()\n",
- "p_nullspace = VectorSpaceBasis(constant=True) # Constant nullspace for pressure n\n",
- "Z_nullspace = MixedVectorSpaceBasis(Z, [Z.sub(0), p_nullspace]) # Setting mixed nullspace\n",
- "\n",
- "# Generating near_nullspaces for GAMG:\n",
- "nns_x = Function(V).interpolate(Constant([1., 0.]))\n",
- "nns_y = Function(V).interpolate(Constant([0., 1.]))\n",
- "V_near_nullspace = VectorSpaceBasis([nns_x, nns_y, x_rotV])\n",
- "V_near_nullspace.orthonormalize()\n",
- "Z_near_nullspace = MixedVectorSpaceBasis(Z, [V_near_nullspace, Z.sub(1)])\n",
- "\n",
- "# To access functions with their appropriate names\n",
- "u, p = z.split() # Do this first to extract individual velocity and pressure fields.\n",
- "# Next rename for output:\n",
- "u.rename(\"Velocity\")\n",
- "p.rename(\"Pressure\")\n",
- "\n",
- "# Setup problem and solver objects so we can reuse (cache) solver setup\n",
- "stokes_problem = NonlinearVariationalProblem(F_stokes, z) #, bcs=[bcvt]) # velocity BC's handled through Nitsche\n",
- "#stokes_problem = NonlinearVariationalProblem(F_stokes, z, bcs=[bcvt]) # velocity BC's handled through Nitsche\n",
- "stokes_solver = NonlinearVariationalSolver(\n",
- " stokes_problem,\n",
- " solver_parameters=stokes_solver_parameters,\n",
- " appctx={\"mu\": mu},\n",
- " nullspace=Z_nullspace,\n",
- " transpose_nullspace=Z_nullspace,\n",
- " near_nullspace=Z_near_nullspace\n",
- ")\n",
- "energy_problem = NonlinearVariationalProblem(F_energy, Tnew, bcs=[bctb, bctt])\n",
- "energy_solver = NonlinearVariationalSolver(energy_problem, solver_parameters=energy_solver_parameters)\n"
+ "energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)\n",
+ "stokes_solver = StokesSolver(z, T, approximation, bcs=stokes_bcs,\n",
+ " cartesian=False,\n",
+ " nullspace=Z_nullspace, transpose_nullspace=Z_nullspace,\n",
+ " near_nullspace=Z_near_nullspace)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Main time stepping\n",
- "We can now initiate the time-loop, with the Stokes system solved seperately. These `solve` calls once again convert symbolic mathematics into computation. In the time loop, set here to run until we reach a steady-state (see `maxchange`), we compute the RMS velocity and surface Nusselt number for diagnostic purposes, and print these results every 50 timesteps. "
+ "We can now initiate the time-loop, with the Stokes and energy systems solved seperately. These `solve` calls once again convert symbolic mathematics into computation. In the time loop, set here to run for 2000 time-steps, we compute the RMS velocity and surface Nusselt number (using definitions from Jarvis, 1993) for diagnostic purposes, and print these results every 50 timesteps."
]
},
{
@@ -454,18 +293,14 @@
"metadata": {},
"outputs": [],
"source": [
- "# Now perform the time loop:\n",
- "for timestep in range(0, max_timesteps):\n",
+ "no_timesteps = 2000\n",
+ "time = 0\n",
"\n",
- " current_delta_t = delta_t\n",
- " if timestep != 0:\n",
- " delta_t.assign(compute_timestep(u, current_delta_t)) # Compute adaptive time-step\n",
- " time += float(delta_t)\n",
+ "for timestep in range(0, no_timesteps):\n",
+ " dt = t_adapt.update_timestep(u)\n",
+ " time += dt\n",
"\n",
- " # Solve Stokes sytem:\n",
" stokes_solver.solve()\n",
- "\n",
- " # Temperature system:\n",
" energy_solver.solve()\n",
"\n",
" # Compute diagnostics:\n",
@@ -473,45 +308,45 @@
" f_ratio = rmin/rmax\n",
" top_scaling = -1.3290170684486309 # log(f_ratio) / (1.- f_ratio)\n",
" bot_scaling = -0.7303607313096079 # (f_ratio * log(f_ratio)) / (1.- f_ratio)\n",
- " nusselt_number_top = (assemble(dot(grad(Tnew), n) * ds(top_id)) / assemble(Constant(1.0, domain=mesh)*ds(top_id))) * top_scaling\n",
- " nusselt_number_base = (assemble(dot(grad(Tnew), n) * ds(bottom_id)) / assemble(Constant(1.0, domain=mesh)*ds(bottom_id))) * bot_scaling\n",
- " energy_conservation = abs(abs(nusselt_number_top) - abs(nusselt_number_base))\n",
- " average_temperature = assemble(Tnew * dx) / domain_volume\n",
- "\n",
- " # Calculate L2-norm of change in temperature:\n",
- " maxchange = sqrt(assemble((Tnew - Told)**2 * dx))\n",
- "\n",
- " if timestep % 100 == 0 or maxchange < steady_state_tolerance:\n",
- " PETSc.Sys.Print(f\"u_rms={u_rms}, Nu_t={nusselt_number_top}, Nu_b={nusselt_number_base}, maxchange={maxchange}\")\n",
+ " nu_top = (assemble(dot(grad(T), n) * ds_t) / assemble(1 * ds_t(domain=mesh))) * top_scaling\n",
+ " nu_base = (assemble(dot(grad(T), n) * ds_b) / assemble(1 * ds_b(domain=mesh))) * bot_scaling\n",
"\n",
- " # Leave if steady-state has been achieved:\n",
- " if maxchange < steady_state_tolerance:\n",
- " log(\"Steady-state achieved -- exiting time-step loop\")\n",
- " break\n",
- "\n",
- " # Set Told = Tnew - assign the values of Tnew to Told\n",
- " Told.assign(Tnew)"
+ " if timestep % 50 == 0:\n",
+ " print(timestep, dt, u_rms, nu_top, nu_base)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can now visualise the final temperature field:"
]
},
{
"cell_type": "code",
"execution_count": null,
- "metadata": {},
+ "metadata": {
+ "tags": [
+ "nbval-ignore-output"
+ ]
+ },
"outputs": [],
"source": [
- "# Let's visualise the initial temperature field now\n",
- "# NBVAL_IGNORE_OUTPUT\n",
- "fig, axes = plt.subplots()\n",
- "collection = tripcolor(Tnew, axes=axes, cmap='coolwarm')\n",
- "fig.colorbar(collection);"
+ "File(\"temp.pvd\").write(T)\n",
+ "temp_data = pv.read(\"temp_0.vtu\")\n",
+ "plotter = pv.Plotter(notebook=True)\n",
+ "plotter.add_mesh(temp_data)\n",
+ "plotter.camera_position = \"xy\"\n",
+ "\n",
+ "plotter.show(jupyter_backend=\"static\", interactive=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "### Exercise 7.2\n",
- "- Change the top boundary conditions from free-slip to no-slip (zero velocity) in the following script. The appropriate boundary condition can be set by using `DirichletBC`. Notice the change in `nullspaces` for this problem. "
+ "### Exercise\n",
+ "2. Change the top boundary condition from free-slip to zero-slip (zero velocity). Note that with a zero-slip boundary condition on either boundary, there is no velocity `nullspace` for this problem. How do your results change, is this what you would expect?"
]
}
],
@@ -519,11 +354,6 @@
"interpreter": {
"hash": "2834c7709b7175be062475c9c8d239d84f48b92a4db284519cf4bc00989e411b"
},
- "kernelspec": {
- "display_name": "Python 3.9.12 ('firedrake')",
- "language": "python",
- "name": "python3"
- },
"language_info": {
"codemirror_mode": {
"name": "ipython",
@@ -534,10 +364,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.9.12"
- },
- "orig_nbformat": 4
+ "version": "3.11.4"
+ }
},
"nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
}
diff --git a/08-GD-2D-pde-constrained1.ipynb b/08-GD-2D-pde-constrained1.ipynb
new file mode 100644
index 0000000..5dfa6aa
--- /dev/null
+++ b/08-GD-2D-pde-constrained1.ipynb
@@ -0,0 +1,619 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "f22c29d2-19b8-440d-bfdf-df258322a99b",
+ "metadata": {},
+ "source": [
+ "# PDE Constrained Optimisation with G-ADOPT\n",
+ "## Example: Inversion for Initial Condition with given Final State\n",
+ "\n",
+ "In this tutorial example we will see how we can use the PDE constrained optimisation functionality of G-ADOPT/Firedrake to optimize one of the inputs to a PDE for a specified desired outcome. We will use a time-dependent advection-diffussion equation as our PDE and see how, for a given final state of the solution we can retrieve what the initial condition should be."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "742b8007-fa95-4ade-8e9c-7b9759052357",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Load Firedrake on Colab\n",
+ "try:\n",
+ " import firedrake\n",
+ "except ImportError:\n",
+ " !wget \"https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh\" -O \"/tmp/firedrake-install.sh\" && bash \"/tmp/firedrake-install.sh\"\n",
+ " import firedrake"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "04c1573f-624c-40f0-b093-170d6a81581d",
+ "metadata": {},
+ "source": [
+ "The usual imports, except we also import the inverse module of gadopt `gadopt.inverse`, which provides\n",
+ "the interface to the inversion functionality we need, but also switches\n",
+ "on the \"tape\" which automatically registers all operations that form\n",
+ "part of the forward model, which is used to automatically form the adjoint (backward) model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "2a4ba5bf-69c1-4649-ae30-a9e4d4f96cf1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Usual imports\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "from gadopt import *\n",
+ "# Starts the tape\n",
+ "from gadopt.inverse import *"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1e5717ab-e34a-445c-bb5f-91f08024f488",
+ "metadata": {},
+ "source": [
+ "## Create Synthetic Twin Experiment With Final State for a Known Initial Condition\n",
+ "We first run a advection-diffusion model with known initial condition. Of that model we only store the solution at the final timestep, which we will use in our inversion experiment later on as the target final state.\n",
+ "\n",
+ "Setup of the model, with mesh, functionspaces, prescribed velocity field, boundary conditions, etc. We use the `EnergySolver` of G-ADOPT to set up an energy equation in the Boussinesq Approximation, which is just a scalar advection-diffusion equation for temperature."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "093bd2b5-d76f-48be-aebb-fe042fac56a3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mesh = UnitSquareMesh(20, 20)\n",
+ "left, right, bottom, top = 1, 2, 3, 4 # Boundary IDs\n",
+ "\n",
+ "V = VectorFunctionSpace(mesh, \"CG\", 2) # function space for velocity\n",
+ "Q = FunctionSpace(mesh, \"CG\", 1) # function space for the scalar (Temperature)\n",
+ "T = Function(Q, name='Temperature')\n",
+ "\n",
+ "# set up velocity field: anti-clockwise rotation around (0.5, 0.5)\n",
+ "x, y = SpatialCoordinate(mesh)\n",
+ "u = interpolate(as_vector((-y+0.5, x-0.5)), V)\n",
+ "u.rename('Velocity')\n",
+ "\n",
+ "# Rayleigh number Ra is not actually used here, kappa is the diffusivity\n",
+ "approximation = BoussinesqApproximation(Ra=1, kappa=1e-2)\n",
+ "temp_bcs = {} # all closed boundaries by default\n",
+ "delta_t = .1 # timestep\n",
+ "energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "90124c96-15f4-47ec-8d1f-283d20ff24af",
+ "metadata": {},
+ "source": [
+ "Setup the initial condition: a Gaussian at $(0.75, 0.5)$. This will be the initial condition we will try to invert for later on."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "51499bd4-2523-4128-aa7f-a9dde4ea7b7e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x0, y0 = 0.75, 0.5\n",
+ "w = .1\n",
+ "r2 = (x-x0)**2 + (y-y0)**2\n",
+ "T.interpolate(exp(-r2/w**2));"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "214a7207-8569-4d52-a075-9803512fc030",
+ "metadata": {},
+ "source": [
+ "A little plotting helper function to plot the Temperature with specified upper bound for the colour scale."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "6f0355a8-a2ee-4fa5-8f30-1ac4dc1f2d3f",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAAGzCAYAAAA/oi4aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAv20lEQVR4nO3de3RV5YH//08unBNCSIgFkkADERlFRQFBMmgZSk3Foji0dRUvhcBSGCu6KvFGREkQJYwXZI2gVAShVgXrqMtvYbAaZayF/qhAOloFBgmCaAJUSDBADkn27w8nx5xc9z5nn8uTvF9rZa1ms/fZT3aRd55nn0ucZVmWAACAceKjPQAAABAcIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIo5OIS4uTsXFxbb2zcnJ0fTp0x2fY//+/YqLi9OaNWscHxsLNm/erLi4OG3evNm/bfr06crJybF1fHFxseLi4sIzOABBIeKICWvWrFFcXJw+/PBDVx5vy5YtKi4u1vHjx115vGBUVlbq7rvv1pAhQ5ScnKwePXpo5MiRevjhh6M6rvacPHlSxcXFAaEHELsSoz0AwA2nTp1SYuJ3f523bNmiBQsWaPr06erVq1fAvrt371Z8fHh/f/3rX/+qiRMn6ptvvtEvf/lLjRw5UpL04YcfavHixXr//ff1xz/+MaxjsGPlypVqaGjwf3/y5EktWLBAkvTDH/4wYN8HHnhAc+fOjeTwAHSAiKNTSEpKsr2v1+sN40ik48eP66c//akSEhK0c+dODRkyJODPH3nkEa1cuTKsY7CrW7dutvdNTEwM+EUJQPSxnI6YNX36dKWkpOjQoUOaPHmyUlJS1KdPH919992qr68P2LfpPfHi4mLdc889kqSzzz5bcXFxiouL0/79+yW1vCf+9ddf6+6779ZFF12klJQUpaam6ic/+Yn+9re/BTXu3/zmNzp06JCWLFnSIuCSlJGRoQceeCBg29NPP60LL7xQXq9X/fr10+zZs1ssuf/whz/U0KFD9cknn2j8+PFKTk5W//799eijj7Y4xxdffKHJkyerR48e6tu3r+bMmaPa2toW+zW9J75//3716dNHkrRgwQL/dWt6XZvfE6+rq9PChQt1zjnnyOv1KicnR/fff3+Lc+Xk5Oiaa67RBx98oNGjRyspKUmDBg3Sb3/723avJYD2EXHEtPr6ek2YMEHf+9739Pjjj2vcuHF64okn9Oyzz7Z5zM9+9jPdcMMNkqQnn3xSL7zwgl544QV/oJrbt2+f3njjDV1zzTVasmSJ7rnnHn300UcaN26cvvzyS8djfvPNN9W9e3ddd911tvYvLi7W7Nmz1a9fPz3xxBP6+c9/rt/85je68sordebMmYB9jx07pquuukrDhg3TE088oSFDhui+++7Tf/3Xf/n3OXXqlK644gq99dZbuv322zVv3jz96U9/0r333tvuOPr06aNnnnlGkvTTn/7Uf91+9rOftXnMLbfcovnz5+uSSy7Rk08+qXHjxqmkpETXX399i3337t2r6667Tj/+8Y/1xBNPKD09XdOnT9ff//53W9cJQCssIAY8//zzliTrr3/9q39bfn6+Jcl66KGHAvYdMWKENXLkyIBtkqyioiL/94899pglySovL29xroEDB1r5+fn+70+fPm3V19cH7FNeXm55vd6Ac5eXl1uSrOeff77dnyU9Pd0aNmxYu/s0Onz4sOXxeKwrr7wyYAzLli2zJFmrV6/2bxs3bpwlyfrtb3/r31ZbW2tlZmZaP//5z/3bli5dakmyXnnlFf+2mpoaa/DgwZYk67333vNvz8/PtwYOHOj//siRIy2uZaOioiKr6T8ZZWVlliTrlltuCdjv7rvvtiRZ7777rn/bwIEDLUnW+++/H/Cze71e66677urgKgFoCzNxxLxbb7014PuxY8dq3759rj2+1+v1P9Gtvr5e//jHP5SSkqLzzjtPO3bscPx41dXV6tmzp61933nnHfl8Pt15550BT7abOXOmUlNTtWHDhoD9U1JS9Mtf/tL/vcfj0ejRowOux8aNG5WVlRWwEpCcnKxZs2Y5/lnas3HjRklSQUFBwPa77rpLklqM/YILLtDYsWP93/fp00fnnXeeq/9fAl0NEUdMS0pKarEMnp6ermPHjrl2joaGBj355JP6p3/6J3m9XvXu3Vt9+vTR//zP/6iqqsrx46WmpurEiRO29v38888lSeedd17Ado/Ho0GDBvn/vNH3v//9Fvelm1+Pzz//XIMHD26xX/NzhOrzzz9XfHy8Bg8eHLA9MzNTvXr1ajH2AQMGtHgMt/+/BLoaIo6YlpCQEPZzLFq0SAUFBfqXf/kX/e53v9Nbb72lt99+WxdeeGHAy6/sGjJkiPbs2SOfz+f6WNu6HpZluX4uu+y+AUwsjh0wHRFHp+TkncVeffVVjR8/XqtWrdL111+vK6+8Unl5eUG/IcukSZN06tQp/ed//meH+w4cOFDSt69db8rn86m8vNz/504MHDhQn332WYs4Nj9Ha5xct4EDB6qhoUH/+7//G7C9srJSx48fD2rsAJwh4uiUevToIUm2QpyQkNAieL///e916NChoM596623KisrS3fddZf27NnT4s8PHz6shx9+WJKUl5cnj8ej//iP/wgYw6pVq1RVVaWrr77a8fknTpyoL7/8Uq+++qp/28mTJ9t9Rn+j5ORkSfau28SJEyVJS5cuDdi+ZMkSSQpq7ACc4Z0b0Ck1vkPavHnzdP3116tbt26aNGmSP+5NXXPNNXrooYc0Y8YMXXbZZfroo4/04osvatCgQUGdOz09Xa+//romTpyo4cOHB7xj244dO/Tyyy9rzJgxkr59cldhYaEWLFigq666Stdee612796tp59+WpdeemnAk9jsmjlzppYtW6Zp06Zp+/btysrK0gsvvOAPdHu6d++uCy64QOvXr9e5556rs846S0OHDtXQoUNb7Dts2DDl5+fr2Wef1fHjxzVu3Dht27ZNa9eu1eTJkzV+/HjHYwfgDBFHp3TppZdq4cKFWrFihTZt2qSGhgaVl5e3GvH7779fNTU1eumll7R+/Xpdcskl2rBhQ0hvMZqbm6uPP/5Yjz32mDZs2KAXXnhB8fHxOv/88zV37lzdfvvt/n2Li4vVp08fLVu2THPmzNFZZ52lWbNmadGiRY7eUa1RcnKySktLdccdd+ipp55ScnKybrrpJv3kJz/RVVdd1eHxzz33nO644w7NmTNHPp9PRUVFrUa8cd9BgwZpzZo1ev3115WZmanCwkIVFRU5HjcA5+IsnlUCAICRuCcOAIChiDgAAIYi4gAAGMpxxN9//31NmjRJ/fr1U1xcnN54440Oj9m8ebMuueQSeb1eDR48WGvWrAliqAAAxKZotdFxxGtqajRs2DAtX77c1v7l5eW6+uqrNX78eJWVlenOO+/ULbfcorfeesvxYAEAiEXRamNIz06Pi4vT66+/rsmTJ7e5z3333acNGzbo448/9m+7/vrrdfz4cW3atCnYUwMAEJMi2cawv05869atysvLC9g2YcIE3XnnnW0eU1tbq9raWv/3DQ0N+vrrr/W9733P0dtCAgCiz7IsnThxQv369Qv4tL5IOn36tOPPM7Asq0VzvF6vvF5vyOMJpo2tCXvEKyoqlJGREbAtIyND1dXVOnXqlLp3797imJKSEi1YsCDcQwMARNDBgwf1/e9/P+LnPX36tAYM7KEjh519oFFKSoq++eabgG1FRUUqLi4OeUzBtLE1MfmObYWFhQGfUVxVVaUBAwZoXK8blBjnieLIAABO1Vk+/ffxl9WzZ8+onN/n8+nI4QZt2dZHKSn2VnO/+cbSZaOP6ODBg0pNTfVvd2MW7qawRzwzM1OVlZUB2yorK5WamtrmbxptLVckxnmUGE/EAcAo/zcBjvbt0JSUOPXsaXc5/9tBp6amBkTcLcG0sTVhvzkxZswYlZaWBmx7++23/R8AAQBAV+NWGx1H/JtvvlFZWZnKysokffs0+bKyMh04cEDSt0vh06ZN8+9/6623at++fbr33nu1a9cuPf3003rllVc0Z84cp6cGACAmRauNjiP+4YcfasSIERoxYoQkqaCgQCNGjND8+fMlSV999ZV/0JJ09tlna8OGDXr77bc1bNgwPfHEE3ruuec0YcIEp6cGACAmRauNRnyKWXV1tdLS0nRFej73xAHAMHUNPpUeW6uqqqqw3F/uSGND/ueTvrbviZ840aCLLzgctTHbxXunAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIYi4gAAGIqIAwBgKCIOAIChEqM9AAAAImHzqUHqnmAve6dO1Uk6HN4BuYCZOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGCivjy5cuVk5OjpKQk5ebmatu2be3uv3TpUp133nnq3r27srOzNWfOHJ0+fTqoAQMAEKsi3UfHEV+/fr0KCgpUVFSkHTt2aNiwYZowYYIOHz7c6v4vvfSS5s6dq6KiIn366adatWqV1q9fr/vvv9/pqQEAiFnR6KPjiC9ZskQzZ87UjBkzdMEFF2jFihVKTk7W6tWrW91/y5Ytuvzyy3XjjTcqJydHV155pW644YYOfzsBAMAk0eijo4j7fD5t375deXl53z1AfLzy8vK0devWVo+57LLLtH37dv+g9u3bp40bN2rixIltnqe2tlbV1dUBXwAARFrzFtXW1ra6X6T62Fyig59FR48eVX19vTIyMgK2Z2RkaNeuXa0ec+ONN+ro0aP6wQ9+IMuyVFdXp1tvvbXd5YKSkhItWLDAydAAAGjXluOD5anz2NrX941P0l+UnZ0dsL2oqEjFxcUt9o9UH5sL+7PTN2/erEWLFunpp5/Wjh079Nprr2nDhg1auHBhm8cUFhaqqqrK/3Xw4MFwDxMAgBYOHjwY0KPCwkLXHjuYPjbnaCbeu3dvJSQkqLKyMmB7ZWWlMjMzWz3mwQcf1NSpU3XLLbdIki666CLV1NRo1qxZmjdvnuLjW/4e4fV65fV6nQwNAADXpaamKjU1tcP9ItXH5hzNxD0ej0aOHKnS0lL/toaGBpWWlmrMmDGtHnPy5MkWA0lISJAkWZbl5PQAAMSkaPXR0UxckgoKCpSfn69Ro0Zp9OjRWrp0qWpqajRjxgxJ0rRp09S/f3+VlJRIkiZNmqQlS5ZoxIgRys3N1d69e/Xggw9q0qRJ/sECAGC6aPTRccSnTJmiI0eOaP78+aqoqNDw4cO1adMm/838AwcOBPxm8cADDyguLk4PPPCADh06pD59+mjSpEl65JFHnJ4aAICYFY0+xlkGrGlXV1crLS1NV6TnKzHe3jMLAQCxoa7Bp9Jja1VVVWXr/rLbGhtyfekv5Umx/+z0dVf8Lmpjtov3TgcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQyVGewAAAETCx//IUsIpr61960/Whnk07mAmDgCAoYg4AACGIuIAABiKiAMAYCgiDgCAoYg4AACGIuIAABiKiAMAYCgiDgCAoYg4AACGIuIAABiKiAMAYCgiDgCAoYg4AACGIuIAABiKzxMH0KVYZ/f3/++48kNRHAkQOiIOoNNrGu62thN0mIiIAzBWW3EO12MResQaIg7ACG4G2+0xEHdECxEHEHNiIdhONI6XmCPSiDiAqDIt2O2xzu5PyBFRRBxAxHWmcDdHyBFJvE4cQER15oA3ss7u3yV+TkQfM3EAXULNgB4ttvU4UBPWc3KvHOFGxAFETDhnp61F2skx4Qw6MUe4EHEAERGOgAcTbruPFY6oE3O4jYgDMI6b8XZyDrfCTszhFiIOIOzcmIVHItzBjiHYuPO2rwgVEQcQVsEEPBaC7UTjeEOZqTM7RzCIOABXhTrrNi3gTdUM6BHykjuvM4cTvE4cgGtiPeDf9EsI+AoHN34GXmMOu5iJA+i0Ogp10z9P+bLetfO6tbzOjBwdCWomvnz5cuXk5CgpKUm5ubnatm1bu/sfP35cs2fPVlZWlrxer84991xt3LgxqAEDiE2xMgsPdqYdjll6qD8TM3LzRLqPjmfi69evV0FBgVasWKHc3FwtXbpUEyZM0O7du9W3b98W+/t8Pv34xz9W37599eqrr6p///76/PPP1atXL6enBtBJuRFwt5fH3Zqlh3qfnBm5OaLRR8cRX7JkiWbOnKkZM2ZIklasWKENGzZo9erVmjt3bov9V69era+//lpbtmxRt27dJEk5OTlOTwsALYTrvnZ75wkm6KEurxNyM0Sjj46W030+n7Zv3668vLzvHiA+Xnl5edq6dWurx7z55psaM2aMZs+erYyMDA0dOlSLFi1SfX3b/yHU1taquro64AtA7Ap22bdmQA/Hs/BwPzHNyfmdjsHkZ953Vc1bVFtb2+p+kepjc45m4kePHlV9fb0yMjICtmdkZGjXrl2tHrNv3z69++67uummm7Rx40bt3btXt912m86cOaOioqJWjykpKdGCBQucDA1AFIQSb6fcjPbJTCm5wp3HahyX3Rl6sMvrzMZDd7wyRfHdk2zt23Dq25lxdnZ2wPaioiIVFxe32D9SfWwu7M9Ob2hoUN++ffXss88qISFBI0eO1KFDh/TYY4+1OcjCwkIVFBT4v6+urm5xIQFEVyTexMWNcJ/M7Hi7G0H/pl9C2EOOyDt48KBSU1P933u9XtceO5g+Nuco4r1791ZCQoIqKysDtldWViozs/X/UrKystStWzclJHz3H+P555+viooK+Xw+eTyeFsd4vV5XLxQA90TiGdPBxrutYDs9LtiohzvkzMYjLzU1NSDibYlUH5tzdE/c4/Fo5MiRKi0t9W9raGhQaWmpxowZ0+oxl19+ufbu3auGhgb/tj179igrK8vWAAHEjlACbmcW7vQ+88nMwC+3hPK4TsbPPfLOI1p9dPw68YKCAq1cuVJr167Vp59+ql/96leqqanxPxtv2rRpKiws9O//q1/9Sl9//bV+/etfa8+ePdqwYYMWLVqk2bNnOz01gCgK5ww82HhHitOohzPkvHY8dkWjj47viU+ZMkVHjhzR/PnzVVFRoeHDh2vTpk3+m/kHDhxQfPx3vxtkZ2frrbfe0pw5c3TxxRerf//++vWvf6377rvP6akBGKqtUDldNg813LWZZ1ps81Z0c/w4du6nc4+864lGH+Msy7Jc/0lcVl1drbS0NF2Rnq/EeJbggWhweynd6czbqdaC3Z5gYt5UWzF38rpyuyE37b54XYNPpcfWqqqqytb9Zbc1NuT7Ty1w8Oz00/rijqKojdkuPgAFQFiFct/XybJ5beaZgC+ngj2uUVvjDMfSOkvqaMQHoAAIi1DjbVco4e3o8ZzOztt6/bmTpXXACWbiADrkdObXUcDbmp0GM/MORnrmibCdI9QZObNxOMFMHIBrYnH23Vawm24/VtHT0fk6mqGHOiPniW6wi5k4AFfYDXgo78LmZGacnnnC9ozbyb5Nx9HeeCIxI2c2DmbiAKKuo1m4GzNvp8d2NDtvqjbzTKuzc+6RI9yYiQNol5uzPaez8HDNvMPxeOGYkTMbR0eYiQMIWTjuhTuJt13D+nzZYtvfjvSz9fh2ZubhmJHbuT/Oe6p3XczEAUSE25//7WSmPKzPl60GvKM/a+18HZ2TGTkiiZk4gJjT3iw81Jl3R/t2NDNvOoa2ZufcI0ekMBMHEHbtvS7cCTdm3naPdTI7b43bM3Jm42gNEQfQJjtRcPvjNNuKn5OAu8XuY7n5hDrACSIOICrc/ihRJ7Pvsb32aGyvPbYfN1huvyUs0BwRBxBxTp+R3tFM10lom8Y7EiFvTbC/wLCkjuaIOICg2YlKqM9KD1fA29sWzHncWFJ3+xn86PyIOICYEMzSs1vL59GakQOhIuIAghLsLNzJUnJ7s1snAXdzv/bO29p4uS+OcCLiAKLOaejcDrjT/d2YkYf6ASmARMQBxKhQ7zE7DbjT41haRywg4gBaFeqznO3OKN2ehTt5+Vh7jxHsWHjNOCKJiANwzO03eGmurRDaCbhb3JyRu3lfPNzXHmYh4gBc59bbrDbl1vJ1XvJe/1dHgv2lwO5snPviCBURBxBTgl2OthPc5uF2K+SRno13hDd86TqIOICY58YyelvBDteMnHvjiAQiDsBVTpaC3ZidhhJwu39u5zxuz8bbu47cF0cjIg4gprUXRzcC7mQ/N58418jtD4JB10LEAbQQjnuqdmLl9hK03YA72d9JyFlSR7gRcQBG6iimTgPu5Li2zh1LS+o8ua1rIOIAjBOOZe2m3JyR25mNs6SOYCVGewAAuiY7M9JgXxse7Cy8+WO8c3JwyI+D2OGp7KaEpG629q0/XR/m0biDmTiAqPBWdPyP6d+O9Avqsd2ILwGHCYg4gIhIruh4n2MVPW091p+On9vhPuGOsJ0xAOFGxAG0EFd+KOhjU760vwwZymw8nCHv6DgnAbfzi4mdX3CcCuX/Q5iDiANwrMeBmmgPQVJ4Qh5KwINd/geCRcQBRIybS+qN3Aw598FhGiIOIKrceIKbGyG3E3A3ZuF2ft5GTm5NoGsi4gCC0t6SeqjxcTobl0ILeagBB6KFiAOIOrdebuY05O+cHOxKwNsaWyhPamMWDjuIOICIsvtM7GBm45L9kNu9/x2rM/D2VkJ4ZnrXQcQBtMpOCMK5pN4au/edIxneUGbhQKiIOICYYPcJX5EK+Z+On+v6Mrrdn5GldNhFxAHErLZiGO6Q2zmO14QjFhBxAGHT1ozSjXcoC1fIQ53BO5mFt3YdmIXDCSIOIGa0Frr27i27GXI7y+dOzxsusfKOeYg+Ig7AaG6E3Mnsu73zhXovHHCKiAMIq0gsD4cScrcC7lSwS+nMwtEUEQcQFU7ui9t5uVYwIXeyfN7R4zMLRzQQcQBtCvW14sEIJXpOQh7N+9/h+OhRdE1EHIAR7L55ipvRtftYkZqF2/mFiXdr61oSoz0AAJ1fypf1+qZfQovtyRXSycyW+3sruqk280yL7Y2xTM880e75msZ3WJ8vHY7W2S8CTgPOe6XDTczEAbTLrSV1p68Zb28We6yip6OZud0oO9nXyRiAcCHiADoUiyGXvgup3Se+tRXoSMWbWTjcxnI6AFviyg/JOrt/u/v0OFCjmgE92t3HraX15oJZanfCSbjD8Yx07oejNczEAUSc2zPyptxe5nbyeN6KbtwLR0QRcQC2ufmSs3CGXHK21N7e8Xa0F28gnFhOB+CIW8vqUnBL603ZWWaX7C+1Ow2+3XCHOgvnXdrQFiIOICzCFfKmnEbdrWX2SM26ncSb++FdE8vpABxzOxhufWRpuJe1nT5+ckXws3Bm37CDiAMIGychcvOzxxtj62bQ3Yp3R3ocqCHgsI3ldABhZXdZXQptab0tjfG1e/+8+XFO2A13W7+wEG84RcQBBMXOE9wauRVyKTwxD3XGzgeaIFqIOICIcCPkUstgOo26m0vsTuPd3n3wUGbhPKmt6+KeOICgOY2HG/fIm2u8/xzJ2bDT86V8Wc+buiAsgor48uXLlZOTo6SkJOXm5mrbtm22jlu3bp3i4uI0efLkYE4LoBMIR8gbNQ26W1Fv/phO7nvbjXews/C48kPMwmNMpPvoOOLr169XQUGBioqKtGPHDg0bNkwTJkzQ4cOH2z1u//79uvvuuzV27FinpwQQw4KJSDhD3pTT+AYb7KYiNesm3rEnGn10HPElS5Zo5syZmjFjhi644AKtWLFCycnJWr16dZvH1NfX66abbtKCBQs0aNCgDs9RW1ur6urqgC8AnUukQt5Ua5F2a+buNN6NLyVzOgtn9h1ZzVtUW1vb5r6R6GNzjiLu8/m0fft25eXlffcA8fHKy8vT1q1b2zzuoYceUt++fXXzzTfbOk9JSYnS0tL8X9nZ2U6GCSDCgo2K05DH4n3lYOMdDOIdmuTDHf8i5//6v8lzdnZ2QI9KSkpafexI9bE5R89OP3r0qOrr65WRkRGwPSMjQ7t27Wr1mA8++ECrVq1SWVmZ7fMUFhaqoKDA/311dTUhBzopJ89al1rOytt6Fnu4hPKLBM9AN8/BgweVmprq/97r9ba6X6T62FxYX2J24sQJTZ06VStXrlTv3r1tH+f1etu8UABik5PXjTfXGDcnMW/UNKrhDHqoqwC8kYuZUlNTAyLulmD72JyjiPfu3VsJCQmqrKwM2F5ZWanMzJYv1vzss8+0f/9+TZo0yb+toaHh2xMnJmr37t0655xzghk3gBgUSsgl57Py5twOeqws3zMLj33R6qOje+Iej0cjR45UaWlpwElLS0s1ZsyYFvsPGTJEH330kcrKyvxf1157rcaPH6+ysjKWyIFOKNTguDVjbfoSL6cxdvv+O7Pwzi9afXS8nF5QUKD8/HyNGjVKo0eP1tKlS1VTU6MZM2ZIkqZNm6b+/furpKRESUlJGjp0aMDxvXr1kqQW2wF0Hm7MyKXgltfb0t4sPZwzbgLedUSjj44jPmXKFB05ckTz589XRUWFhg8frk2bNvlv5h84cEDx8bwRHNDVhRpyKfTl9bbEyjJ5R1hGN0s0+hhnWZbl6iOGQXV1tdLS0nRFer4S4z3RHg4AB0INueTujDySgpmFd8Zw1zX4VHpsraqqqsLyJLGONDZk6KxFSvAk2Tqm3ndaHz97f9TGbBcfgAIgrBqjFGvL625zY9m8MwYc4UXEAURE00CF+lK01kQj8G7e7ybgCAYRBxBxbszOm2se1HBFnSeqIZYQcQBRE46YN3I76uGMN7NwBIuIA4i6cMa8UdMINw96NGfXBByhIOIAYkYkYi6xJI7Ogxd0A4g5XeXjNrvCz4jwYiYOIGa58YYxsYZww01EHEBMMznkBBvhRsQBxLxI3SsPBcFGNBBxAMZoK5SRjjvBRqwg4gCMZyeqHYWeMMNERBxAl9Da274SbpiOiAPocog3OgteJw4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYiogDAGAoIg4AgKGIOAAAhiLiAAAYis8TBwB0CT2+qldit3pb+9adsbdftDETBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEFFfPny5crJyVFSUpJyc3O1bdu2NvdduXKlxo4dq/T0dKWnpysvL6/d/QEAMFWk++g44uvXr1dBQYGKioq0Y8cODRs2TBMmTNDhw4db3X/z5s264YYb9N5772nr1q3Kzs7WlVdeqUOHDjk9NQAAMSsafYyzLMtyMsjc3FxdeumlWrZsmSSpoaFB2dnZuuOOOzR37twOj6+vr1d6erqWLVumadOmtbpPbW2tamtr/d9XV1crOztbV6TnKzHe42S4AIAoq2vwqfTYWlVVVSk1NTXi56+urlZaWppyJy1UYrckW8fUnTmt/+//PaiDBw8GjNnr9crr9bZ6TCT62JyjmbjP59P27duVl5f33QPExysvL09bt2619RgnT57UmTNndNZZZ7W5T0lJidLS0vxf2dnZToYJAEALPb6oUY8DNr++qJEkZWdnB/SopKSk1ceOVB+bS7S9p6SjR4+qvr5eGRkZAdszMjK0a9cuW49x3333qV+/fgE/aHOFhYUqKCjwf984EwcAIJJam4m3JlJ9bM5RxEO1ePFirVu3Tps3b1ZSUttLGu0tVwAAECmpqakRuQVgt4/NOYp47969lZCQoMrKyoDtlZWVyszMbPfYxx9/XIsXL9Y777yjiy++2MlpAQCIadHqo6N74h6PRyNHjlRpaal/W0NDg0pLSzVmzJg2j3v00Ue1cOFCbdq0SaNGjXI0QAAAYl20+uh4Ob2goED5+fkaNWqURo8eraVLl6qmpkYzZsyQJE2bNk39+/f33/z/93//d82fP18vvfSScnJyVFFRIUlKSUlRSkqK4wEDABCLotFHxxGfMmWKjhw5ovnz56uiokLDhw/Xpk2b/DfzDxw4oPj47yb4zzzzjHw+n6677rqAxykqKlJxcbHT0wMAEJOi0UfHrxOPhsbX+PE6cQAwT6y8TvxHI+YqMcHm68TrT+vdnYujNma7eO90AAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMRcQBADAUEQcAwFBEHAAAQxFxAAAMlRjtAQAAEAlx+79SXLzH3r4NvjCPxh3MxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAxFxAEAMBQRBwDAUEQcAABDEXEAAAwVVMSXL1+unJwcJSUlKTc3V9u2bWt3/9///vcaMmSIkpKSdNFFF2njxo1BDRYAgFgW6T46jvj69etVUFCgoqIi7dixQ8OGDdOECRN0+PDhVvffsmWLbrjhBt18883auXOnJk+erMmTJ+vjjz92emoAAGJWNPoYZ1mW5WSQubm5uvTSS7Vs2TJJUkNDg7Kzs3XHHXdo7ty5LfafMmWKampq9Ic//MG/7Z//+Z81fPhwrVixwtY5q6urlZaWpivS85UY73EyXABAlNU1+FR6bK2qqqqUmpoa8fMH05BgxhyNPiba2uv/+Hw+bd++XYWFhf5t8fHxysvL09atW1s9ZuvWrSooKAjYNmHCBL3xxhttnqe2tla1tbX+76uqqiRJdZZPanAyYgBAtNVZPkmSwzljeMZhsyGNY66urg7Y7vV65fV6W+wfqT425yjiR48eVX19vTIyMgK2Z2RkaNeuXa0eU1FR0er+FRUVbZ6npKRECxYsaLH9v4+/7GS4AIAY8o9//ENpaWkRP6/H41FmZqb+u8JZQ1JSUpSdnR2wraioSMXFxS32jVQfm3MU8UgpLCwM+O3k+PHjGjhwoA4cOBCVvwCmqK6uVnZ2tg4ePBiVJStTcJ06xjWyh+tkT1VVlQYMGKCzzjorKudPSkpSeXm5fD6fo+Msy1JcXFzAttZm4dHkKOK9e/dWQkKCKisrA7ZXVlYqMzOz1WMyMzMd7S+1vVyRlpbGfyg2pKamcp1s4Dp1jGtkD9fJnvj46L2qOSkpSUlJSWF7/Ej1sTlHV9Tj8WjkyJEqLS31b2toaFBpaanGjBnT6jFjxowJ2F+S3n777Tb3BwDANFHro+XQunXrLK/Xa61Zs8b65JNPrFmzZlm9evWyKioqLMuyrKlTp1pz58717//nP//ZSkxMtB5//HHr008/tYqKiqxu3bpZH330ke1zVlVVWZKsqqoqp8PtUrhO9nCdOsY1sofrZE9XuU7R6KPjiFuWZT311FPWgAEDLI/HY40ePdr6y1/+4v+zcePGWfn5+QH7v/LKK9a5555reTwe68ILL7Q2bNjg6HynT5+2ioqKrNOnTwcz3C6D62QP16ljXCN7uE72dKXrFOk+On6dOAAAiA28dzoAAIYi4gAAGIqIAwBgKCIOAIChiDgAAIaKmYjzGeX2OLlOK1eu1NixY5Wenq709HTl5eV1eF07A6d/lxqtW7dOcXFxmjx5cngHGCOcXqfjx49r9uzZysrKktfr1bnnntsl/rtzep2WLl2q8847T927d1d2drbmzJmj06dPR2i00fH+++9r0qRJ6tevn+Li4mx9gMfmzZt1ySWXyOv1avDgwVqzZk3Yx9kphf6quNCtW7fO8ng81urVq62///3v1syZM61evXpZlZWVre7/5z//2UpISLAeffRR65NPPrEeeOABxy+QN5HT63TjjTday5cvt3bu3Gl9+umn1vTp0620tDTriy++iPDII8fpNWpUXl5u9e/f3xo7dqz1r//6r5EZbBQ5vU61tbXWqFGjrIkTJ1offPCBVV5ebm3evNkqKyuL8Mgjy+l1evHFFy2v12u9+OKLVnl5ufXWW29ZWVlZ1pw5cyI88sjauHGjNW/ePOu1116zJFmvv/56u/vv27fPSk5OtgoKCqxPPvnEeuqpp6yEhARr06ZNkRlwJxITER89erQ1e/Zs//f19fVWv379rJKSklb3/8UvfmFdffXVAdtyc3Otf/u3fwvrOKPN6XVqrq6uzurZs6e1du3acA0x6oK5RnV1ddZll11mPffcc1Z+fn6XiLjT6/TMM89YgwYNsnw+X6SGGBOcXqfZs2dbP/rRjwK2FRQUWJdffnlYxxlL7ET83nvvtS688MKAbVOmTLEmTJgQxpF1TlFfTm/8DNa8vDz/Njufwdp0f+nbz2Bta//OIJjr1NzJkyd15syZqH2SULgFe40eeugh9e3bVzfffHMkhhl1wVynN998U2PGjNHs2bOVkZGhoUOHatGiRaqvr4/UsCMumOt02WWXafv27f4l93379mnjxo2aOHFiRMZsiq74b3i4RP2jSKP1GaymCeY6NXffffepX79+Lf7j6SyCuUYffPCBVq1apbKysgiMMDYEc5327dund999VzfddJM2btyovXv36rbbbtOZM2dUVFQUiWFHXDDX6cYbb9TRo0f1gx/8QJZlqa6uTrfeeqvuv//+SAzZGG39G15dXa1Tp06pe/fuURqZeaI+E0dkLF68WOvWrdPrr78e1o/jM8mJEyc0depUrVy5Ur179472cGJaQ0OD+vbtq2effVYjR47UlClTNG/ePK1YsSLaQ4spmzdv1qJFi/T0009rx44deu2117RhwwYtXLgw2kNDJxX1mXi0PoPVNMFcp0aPP/64Fi9erHfeeUcXX3xxOIcZVU6v0Weffab9+/dr0qRJ/m0NDQ2SpMTERO3evVvnnHNOeAcdBcH8XcrKylK3bt2UkJDg33b++eeroqJCPp9PHo8nrGOOhmCu04MPPqipU6fqlltukSRddNFFqqmp0axZszRv3ryofp52LGnr3/DU1FRm4Q5F/W8Un1FuTzDXSZIeffRRLVy4UJs2bdKoUaMiMdSocXqNhgwZoo8++khlZWX+r2uvvVbjx49XWVmZsrOzIzn8iAnm79Lll1+uvXv3+n/JkaQ9e/YoKyurUwZcCu46nTx5skWoG3/xsfisKb+u+G942ET7mXWWFZ3PYDWR0+u0ePFiy+PxWK+++qr11Vdf+b9OnDgRrR8h7Jxeo+a6yrPTnV6nAwcOWD179rRuv/12a/fu3dYf/vAHq2/fvtbDDz8crR8hIpxep6KiIqtnz57Wyy+/bO3bt8/64x//aJ1zzjnWL37xi2j9CBFx4sQJa+fOndbOnTstSdaSJUusnTt3Wp9//rllWZY1d+5ca+rUqf79G19ids8991iffvqptXz5cl5iFqSYiLhlRf4zWE3l5DoNHDjQktTiq6ioKPIDjyCnf5ea6ioRtyzn12nLli1Wbm6u5fV6rUGDBlmPPPKIVVdXF+FRR56T63TmzBmruLjYOuecc6ykpCQrOzvbuu2226xjx45FfuAR9N5777X6b03jtcnPz7fGjRvX4pjhw4dbHo/HGjRokPX8889HfNydAZ8nDgCAoaJ+TxwAAASHiAMAYCgiDgCAoYg4AACGIuIAABiKiAMAYCgiDgCAoYg4AACGIuIAABiKiAMAYCgiDgCAof5/gM5Lau6JvNMAAAAASUVORK5CYII=",
+ "text/plain": [
+ "