{ "cells": [ { "cell_type": "markdown", "id": "ee9c5bb7", "metadata": {}, "source": [ ".. index::\n", " pair: Methods; Discontinuous Galerkin\n", "\n", "# Discontinuous Galerkin Methods\n", "\n", "## Advection-Diffusion: Discontinuous Galerkin Method with Upwinding\n", ".. index:: Equations; Advection-Diffusion\n", "\n", ".. index::\n", " pair: Spaces; DG Legendre\n", "\n", "So far we have been using Lagrange spaces of different order to solve our\n", "PDE. In the following we show how to use Discontinuous Galerkin method to\n", "solve an advection dominated advection-diffusion problem:\n", "\\begin{align*}\n", "-\\varepsilon\\triangle u + b\\cdot\\nabla u &= f\n", "\\end{align*}\n", "with Dirichlet boundary conditions. Here $\\varepsilon$ is a small\n", "constant and $b$ a given vector." ] }, { "cell_type": "code", "execution_count": 1, "id": "4e73cc3f", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:28:56.746219Z", "iopub.status.busy": "2024-02-29T12:28:56.745863Z", "iopub.status.idle": "2024-02-29T12:32:00.893318Z", "shell.execute_reply": "2024-02-29T12:32:00.891609Z" } }, "outputs": [], "source": [ "import numpy, math, sys\n", "from matplotlib import pyplot\n", "from dune.grid import structuredGrid as leafGridView\n", "from dune.fem.space import dglegendre as dgSpace\n", "from dune.fem.scheme import galerkin as solutionScheme\n", "from dune.ufl import Constant\n", "from ufl import ( TestFunction, TrialFunction, SpatialCoordinate, triangle, FacetNormal,\n", " dx, ds, grad, div, grad, dot, inner, sqrt, exp, conditional,\n", " as_vector, avg, jump, dS, CellVolume, FacetArea, atan )\n", "\n", "# overlap=1 is needed for parallel computations\n", "gridView = leafGridView([-1, -1], [1, 1], [20, 20], overlap=1)\n", "space = dgSpace(gridView, order=2)\n", "\n", "\n", "# check if computation is parallel\n", "parallel = gridView.comm.size > 1" ] }, { "cell_type": "markdown", "id": "e760a94d", "metadata": {}, "source": [ "\n", ".. todo:: add some more details on the forms for upwinding\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "3409c71d", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:32:00.899723Z", "iopub.status.busy": "2024-02-29T12:32:00.898651Z", "iopub.status.idle": "2024-02-29T12:32:00.925904Z", "shell.execute_reply": "2024-02-29T12:32:00.924443Z" } }, "outputs": [], "source": [ "u = TrialFunction(space)\n", "v = TestFunction(space)\n", "n = FacetNormal(space)\n", "he = avg( CellVolume(space) ) / FacetArea(space)\n", "hbnd = CellVolume(space) / FacetArea(space)\n", "x = SpatialCoordinate(space)\n", "\n", "# diffusion factor\n", "eps = Constant(0.1,\"eps\")\n", "# transport direction and upwind flux\n", "b = as_vector([1,0])\n", "hatb = (dot(b, n) + abs(dot(b, n)))/2.0\n", "# boundary values (for left/right boundary)\n", "dD = conditional((1+x[0])*(1-x[0])<1e-10,1,0)\n", "g = conditional(x[0]<0,atan(10*x[1]),0)\n", "# penalty parameter (take 1 for p=0)\n", "beta = 10 * space.order**2 if space.order > 0 else 1\n", "\n", "aInternal = dot(eps*grad(u) - b*u, grad(v)) * dx\n", "diffSkeleton = eps*beta/he*jump(u)*jump(v)*dS -\\\n", " eps*dot(avg(grad(u)),n('+'))*jump(v)*dS -\\\n", " eps*jump(u)*dot(avg(grad(v)),n('+'))*dS\n", "diffSkeleton += eps*beta/hbnd*(u-g)*v*dD*ds -\\\n", " eps*dot(grad(u),n)*v*dD*ds\n", "advSkeleton = jump(hatb*u)*jump(v)*dS\n", "advSkeleton += ( hatb*u + (dot(b,n)-hatb)*g )*v*dD*ds\n", "form = aInternal + diffSkeleton + advSkeleton" ] }, { "cell_type": "markdown", "id": "0892ddc2", "metadata": {}, "source": [ "We solve the problem as in previous examples" ] }, { "cell_type": "code", "execution_count": 3, "id": "8a5d1598", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:32:00.931065Z", "iopub.status.busy": "2024-02-29T12:32:00.930280Z", "iopub.status.idle": "2024-02-29T12:32:55.058252Z", "shell.execute_reply": "2024-02-29T12:32:55.057135Z" } }, "outputs": [ { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrAS0DASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+ql/frp8cMjwyyJJPHCWj2/IXYIpOSOMsOmTz0q3VTUNNttTgSG6EpRJFlAjmeM7lOVOVIJwQDjpxQBRvfElpYXE8c0NwVhJRpVVSpkEfm+WOc7tnPTHbOasJrdiNNe/u5VsYY3eOU3cip5bKSCCc7e3r0psnh/TJmZpYHkLR+W2+Z2DfJs3HJ5fbxv+9jvWb4ss4LH4c+IoLdWVBpt23zOWJJjckkkkkkk9aBpXdif/AITfwl/0NGif+DCL/wCKo/4Tfwl/0NGif+DCL/4qvkShe/1rl+s+R7v9jL+f8P8Agn13/wAJv4S/6GjRP/BhF/8AFUf8Jv4S/wCho0T/AMGEX/xVfIv8S0rfdNH1nyH/AGKv5/w/4J9c/wDCb+E/+ho0X/wYRf8AxVH/AAm/hP8A6GjRf/BhF/8AFV8kUL3+tL60+xX9hr+f8P8Agn1v/wAJv4T/AOho0X/wYRf/ABVH/Cb+E/8AoaNF/wDBhF/8VXyT/H+FOb7ppfWn2H/YSt8f4f8ABPrX/hN/Cf8A0NGi/wDgfF/8VR/wm3hP/oaNF/8AA+L/AOKr5NoTv9aX1t9i/wDV+P8Az8/D/gn1l/wm3hP/AKGfRf8AwPi/+Ko/4Tbwn/0M+i/+B8X/AMVXycP9Z+FOf7ho+tvsUuHo2b9p+H/BPq//AITbwp/0M+i/+B8X/wAVR/wm3hT/AKGfRf8AwPi/+Kr5UFEf3BU/XH2NFw3G9vafh/wT6r/4Tbwp/wBDPov/AIHxf/FUf8Jr4U/6GfRf/A+L/wCKr5WH+s/CnN/D9aPrr/lKXDMbfxPw/wCCfU//AAmvhT/oZtG/8D4v/iqP+E18Kf8AQz6N/wCB8X/xVfLgoi+4Kn68/wCU0XCsW7e1/D/gn1H/AMJr4U/6GbRv/A+L/wCKpf8AhNfCv/QzaN/4Hxf/ABVfLw/1v/AaeR9z/eFJ49/ylrhONr+1/D/gn09/wmnhX/oZdG/8D4v/AIqj/hM/Cx/5mXRv/A6L/wCKr5sQVNbD92tZvM2vs/iV/qjC/wDF/D/gn0d/wmXhc/8AMyaP/wCB0X/xVL/wmPhj/oY9I/8AA6L/AOKr55g/1z/hV8DmP/fFZSzdx+x+P/AMZ8LRir+1/D/gnvI8W+Gj08Q6Sf8At9j/APiqcvinw84yuvaWR6i8j/xryewHSux8Of8AIMH0P868+txLKn/y7/H/AIB5GNypYaLkp3+R2Vlqmn6l5n2G+tbry8b/ACJlfbnpnB4zg/lU89xDawtNcTRwxL955GCqO3JNc94d/wCQ5qn/AFxt/wCctZnxfGfhfq4/2rf/ANHx19JQre1oqrbdXPGw69s4ra503/CQ6J/0GNP/APAlP8as2uoWd9G72V1BcqhwxgkV8H04NfG3lV9TfDMf8W40Uf8ATFh/4+1OlW9o7WPTx+WvCRUnK9/L/gmnN4ktbbw/Nq80FwiRNKht8KZS8bMrKAGwTlG74wMkgZp1zrwtru5tzp145gtjcl0MWGUEDAy+QTzjIGdrelJ/wi+jtYTWM1o1xbTFmaO5meYBmLFmXex2k725GDzVt9IsXjuI/JKrcRLDJ5bsh2KCAoIIKgZPTHU1seYWYJVngjmUELIoYA9cEZqSora3jtLaO3h3eXGoVd7lzge5JJ/GpaACiioLm8tbJEe6uYYFdgimVwoZj2Gep9qAJ6wfG/8AyIHiP/sF3P8A6KatKfV9NtZ3guNRtIZkQyPHJMqsqAZLEE8ADnPpWT4wnhuvh14hmt5Y5on0u6KvGwZW/dN0IoY47o+SNvuaRV68nrTuf7v60ilufl7+teYfcWVwK/MvJ/OlZflPJ/Oglty/L+tK27afl/WkOysxwT/aP50iJ15PX1pwLf3f1oTdz8vf1pGtldAF+fqenrSsnynk/nQN2/7vb1pzbth+X9aVylFWYoT/AGm/OhE6/M3X1pw3f3P1oTdz8vf1qbm6iroAn7z7zdPWnMnyH5m/OgbvM+729ac+7Yfk/Wpu7miirMUJ/tN+dLHH8g+Zvzpw3f3P1oj3bB8n61N9DdRXMtBBH+8+83T1pzR/d+ZuvrSjd5n3O3rTm3/L8n8Q71N2aKEbPT8xwj/2m/OiKP5B8zfnTxv/ALn60sKuYx+7/Wpb0OhRjzLT8xFizL95vu+tSmH7nzN94d6cqv53+r/h9amYP8n7v+IfxVm5Mq0bP/gjkh/23/OpLeH90Pnf86evmf8APP8A8ep1t5nlL+6/8erCUnYcnG//AA5JBD++k+d+3erwg5j/AHkn3x3qrB5nnSfuvT+Kr4MuY/3X8Y/iFc1Ru5xVZR5Wbtjb9P3sv/fVdd4ets6YP30w4PR/euVsWm4xB/4+K6vQJJxpi4ts9f4x614WKcvxR8tnDXs2bXhlPL1rVRvZv3NucscnrLVL4tjPwz1Uf7dv/wCj46u+GWdta1Uumw+Tb8Zz3lqr8VV3fDjU19ZLf/0fHX6HhNMHH0Pl8tV6tJea/M+bPK9q+m/hsMfDzRx/0yb/ANDavnX7Ka+h/hvNE3gPToo5EkkhVldFYEqd7HB9D9aywE+aTPseKaShRpvzf5HW0VknXohorambO68tJHjki+TehVyjE/NggFT0J4pZ/EFnb6hPZuJS8MJkZ1UFScqNg5yX+dOMfxDvXpnxRq0VW0+9j1LTbW+iR0juYUmVXxuAYAgHBIzz2JqzQAVk67oraxCES5EB8uWFy0e/KSLtbAyMN0wecehzWtRQBl3GkyXkpmuZbeZ/sbWwje3zFl8eYSpblW2qNueADyc1ieJdMfT/AIXa7Zi6IMdhdyFoowq4Ku2wKd21BnaADkADBrr6wfG//Ig+I/8AsF3P/opqT2Kj8SPkbcvrQrLzz3qTFCjr9a8u595yu6GFl3LzSsy7TzT8fMtOYfIaV0WoOzGhl9aEZeee9ShaVF6/WpujZU3dEYZfM69qczLsPNSBf3n4U51+Q1N0aqnKzGhl9aEZfm571OFojX731qeZHQqUroiDL5nXtTnZdh5qYJ+9/wCA050/dmp5lc1VKXKyMMvrSxMuwc1aWEmnwQZjHFZuasbqm+ZFQMvm9f4fSpDs+Xn+IdqurB++/wCA/wBalaL7nH8QqHURXK0nqVU2ev6U+AxiNef0q8IqWCL90vFYuorFOTUtyorR+f1/h9PepXeP5Of4h2qysX+kf8B/rUjxfc/3xUOauZSm7PUhV4/736U+2ePyV+b9KuLF7U61i/crWLmrGc5u+5DBJH58nzenarvmxZi+b+MdjSwR/v5Pw/lV0x8xf74rnqTVziqzfK9TTsriEY+f9DXT6HdwLpagyevY+tY1lH0ro9GGNKT8f5mvGryi2vVHy+by/dMveGJUl1nVWQ5Hk246e8tQ/E0bvAGoD1ltv/R8dWPDn/Ib1T/rjb/zlqP4jjd4GvR6zW3/AKPjr9Fw+mBX+H9D5zLHapSfmvzPDvs59K928ADHgfTV9pB/5EavHvs/tXsfgQY8G2A9DL/6MauDJ5805eh9jxRU56NP1f5F0eGtMGnz2O26NtOxaRTezEklixIO/IyzEnB5709/DujS3JuZdMtZJmTYzyRhi4yp+bP3jlVOTk8Vp0V758WV7GyttNsYLKzhWG2gQRxxr0VR0qxRRQAUUVheJrnU7a2jOm+eGKyEGGHzS0gX92jDBwpPU8YwORmgDdrB8b/8iD4j/wCwXc/+imqlqN7rK6jPJp893NbyWpZLc2BT7OdoO8Oy/O3X5OuT04IqHxDNfy/C7X3lVpJTZXq7rgeU5iAkCsVC/e24OMDPtSexUfiR8vBPc/nQideT19aeA/8As0qK/P3eteTc/Q1BXWg3y/nXlu/enNH8h+Zvzp+196/d709kfYfu1LkbKkrPQaIv9pvzpY4uvzN19amCP/s0saSfN93rWbkdUaKutCMRfvPvN09ae8PyH5m/OplhlMvRfu1M1tL5ZyFqHPXc2VGPK9CAQZ6M3506K2J3fM33vWr6W8von60+GCX5+E+8fWsnVNuWCa0KS2f777zfd9amezxGfmf86uLDL5/RPu+/rUskM3lNwn61k6rutSXKCi9Cutpj+J/zpYLX90vzv+dX1gm9E/Wi3gm8leE/WsXV03JlVV9ioLX9/wDff7vr71I9r9z53+8O9WxBN9o6R/c9/WpJIJv3fEf3x61Dqu+5hKqrMri0/wBt/wA6Le0zEvzv+daCwT+kf60ltDP5C8R9/X1rJ1XbciVVXKi2n+kfff7nr70+S0/1fzyffHerggn+09I/ue/rT5YJ/wB3xH98etS6rvuYyqqzIVtP+mkn50traZhX95J+fvV4Q3HpH+tLaQ3HkLxH39axdV23M5VVcrwWn7+X95J2/i9quGz5i/eyf6wfxU6CG4+0S8Rfw+vpVxobnMXEX+sHrWM6jvuclSqrGjZ2fT97N/31W3plrjS4/wB/MOv8fvWfaxXXHEP61p2IuxpsY/c9/X1rzXJymlfqfOZrP92y/wCGE8vWdVG9m/c2/LHJ6y074gDd4Nuh6z23/o+OmeGPM/tnVfM27vJt/u9OstT+ORu8Jzj1uLb/ANHx1+j09MB/26/yPEy92dN+a/M8w+z16H4Y1OGy8DLLEUnlt3kRolkAwxlYKGPO3qOcdOa5H7Oa7zwnDHL4ShhmjWSN3nV0cZDDzXBBB6ivF4enzVZ+h9LnlTnpw9SG88WPZ2VvMbAO7yTpKiysRGIn2O2Qh+XP8TBVHG4rkVZuddurex1icWUDPp05j2G4IEq+WkgIOw4Y7wNuDz3q+dF0oxrGdMstiPvVfIXCtgDIGOuABn2FPTS9PjknkSwtVeeRZZmWFQZHU5VmOOSCAQT0r6s+aLERkaFGlQJIVBZVbcFPcA8Z+tPoooAKKKKACsLxt/yIXiL/ALBdz/6Kat2sLxr/AMiH4i/7Blz/AOimpPYqHxI+UQV9adGV5571OqVJDETuwP4jXhuSP1SFCTaIMrvXn1qRtuw81eS1/eJn3q21qBEeKwlVSsdHs+VSu/6sUEjQ9T+lTwJH83P8R7Vqpbe1S29tnfx/Ga5pV1YqVSzRmII/P6/w+nvU0gj8o8/pWolr/pPT+D+tTy2v7huPT+dYuuro554h2epnKI/X9KdAI/n+b+I9jW0tr7Utva/6zj+M1zuurGM8Q7rUyV8r7R97+D096kl8ryW+b9D61rra/wCldP4P61LPa4t249P51Drq6OeWIdnqZi+V/e/Q0W3leQvzevY+tbi2vtSWlt/o6cev86xddWM5Yh33MhTF9q+9/B6H1p8pi/d/N/GOxrYW2/0zp/yz/rUk1t/quP8AloKh11dGMsQ7PUzAYf736Gm2rQ+QvzevY+tbwtvamWltm2Tj1/nWTrxsQ8Q77mQGh+1/e/g9D60+Zof3Xzf8tB2NbAtv9N6f8s/60+e3/wBVx/y0FT7aN0Yuu7bmarw/3v0NJZvD9nT5/XsfWtwW3tTbK3/0ZOPX+dZOtHlZnKu77mZBJD9pm+f+72PpVqSWDMPz/wDLQdjV+C3/ANKn4/u/yqxLBjyeP+Wq1lKtHmOedZ2H209uMfP+hqe1vLcWMY8z17H1NXoYgoqvAf8AQ0/H+dZ4W06mx8/mda8GifwvIkusaoyHI8m3HT3lq14zG7w1IPW5tv8A0fHVfw1/yGtU/wCuNv8Azlq54sG7QGHrdW3/AKPjr9KWmXv/AAv8jhwTtGD9Dk/s/tXWeExjw7CPSaf/ANHPWR9nq5o2pLY+Fp5Vid3guJl2srIrMZ2A+YjGMkZIzjmvmuFZ81ap6L8z2cxqc8YnTUVzd54vhsLGzluIoknuLw2pjacKq7ZvKdwxAyASCBjJyOnJDP8AhKrk30tj/ZXl3YkCxrNKyKykTNksU44gb7u4cjnrj7c8k6eioLO6S9sbe7jVlSeNZFDDkBhnn86noAKydd1ptHhV0thOfLlmcNJswka7mwcHLdMDjPqMVrVDc2dreIi3VtDOqMHUSoGCt6jPQ0AZt9qt9b381paadFeSLbNNGiXO1iR91XyuEDHIByeh44OMrxJqsd58MdcurgJau9hdwtG0gO2UK6FM8ZO4Ecde2a3homkh2caXZBnjETN9nTJQAAKeOgAAx6Cs/wAV28Nt8P8AXYLeKOKJdMuQkcahVUeW3QDpSexUPiR8zw2RY/eb86vWthw3zP8AePer8FlPx8qfrV+zsLgh/lj++fWvkauJ03P1mpWhBqyM5dPPmx/O/fvVt9O/ct88nbv71qrYXPnwjbH39fSrkthci3b5Yu3r61xTxWq1POq4yOpmR6af78n51JbabnzPnk++e9b6WFz/AHYv1p1pY3X73iL/AFh9a45YrR6nJUxiuYqab/pWPMk+5/e96nn0zFu37yTt/F71tJY3X23G2L/V+/rU1zY3QtXOIe3r61i8U+ZanLPFqzMpdM/6aS/99UW2mZ8395L/AKw/xV0C2N36Q/rTbWyuj53EP+sPrWDxLs9TCWLVzGXTP9Mx5kv+r/ve9SXOmYtn/eS9v4vetkWd19vxiH/Ve/rT7u0uhauSIe3r61P1l8y1MZYtWZmjTP8AprL/AN9Uyz03Nsn72Xv/ABe9b4tLv0g/Wo7K1uvsicQ9/X1rL6w+XczeLVzHGm/6bjzZf9X/AHvenz6bjyf3s3+sH8Va4trr7f0g/wBV7+tLc290PJ/1H+tX1pe3ldamTxce5RGmf9NZv++qjs9Nzap+9m7/AMXvW35N1/0w/WobKK5+yJ/qO/r61HtZ8u5m8XEzRpv+nY82b/V/3venXGm48n99N/rF/irR2XP9of8ALD/Ve/rS3S3X7n/Uf61fWnz1LozeMiVxpn/Tab/vqo7LTs2qfvpu/wDF71qYuv8Aph+tV7I3QtEwIO/r61H7xxM3jEV4NO/0u4HnTfw/x+1S3FgB5P7+b/Wr/HT4mu/tdxjyP4fX0pLprv8Ac/6j/Wr61ap1HJa/1YwniyybbAwJ5/8AvuqNvb/6In7+bv8Ax+9Wi13/ANMP1qlbG7+yJ/qe/r616OCotM8jF1uaLNPwsmzWNUG5m/c2/LHJ6y1o+J/+QL/29Wv/AKPjrN8LeZ/bGqeZt3eTb/d6dZa0fFHGif8Ab1a/+j46+9avgGv7r/I0wrtSiyptWregRRzaC8MsayRPPcqyOMhgZnyCO4qjurQ8Nf8AIGH/AF83H/o56+Z4Vpclao/JfmdMq3tNC9Fp1lBAYYbO3jiLq5RIlC7hjBwB1G1cH2HpTG0jTHjmjfTrRknfzJVMCkSN/eYY5PuauUV9uSAAAwBgCiiigAooqCa8tbeaGGa5hjlmJESO4DSEddoPXqOlAE9Y3i/nwVrw/wCodcf+i2q8+qafG06vf2qtb484NMo8rPTdzxn3qn4keOfwdq7xurxvp8xVlOQwMZ5B7ipn8LHF2aZ5NbWsPHP/AI6a07C3gxJz/wAtD/Ca69ILaL7qD8aZayIPOwo/1rV+YSlWqLSLPcr8QU5PRmD5FuLqD5v738J9KtXEdsLV/n9P4T61qyzD7Xb4/wBr+VNvZs2j/h/OoWGrSa0PNnnafUiUWo/j/wDHTTbR7QedmQf61v4TV7zveq9tL/ruf+WrUlgKjTvc5pZxfqRia0F//rP+WX90+tOvLq0+yP8AP6fwn1pfN/0/r/yy/rSXsv8Aoj/h/OrWWu6vcxeaX6ln7baj+P8A8dNV7W+th53zj/Wt2NTebUFtJ/rv+urVUct02MnmV+ov2+2/tD74/wBV/dPrRe39v9kf956fwn1pvmf6f/2y/rReSf6I/wCH860WXK60IeYeZZ/tC3/56foar2eoW/2RP3nr2PrU3mVBZyf6In4/zprLlbYh4/zF+32/2/8A1n/LL0PrSXWoW/7n95/y1XsaTzP9P6/8sv60XT/6n/rqtaLAK60J+veZY/tC3/56foags7+3+yJ+89ex9am31BZv/oifj/OmsArbE/XQ+32/2/8A1n/LL+6fWi6v7f8Ac/vP+Wq9jRv/ANP/AO2X9aS5f/U/9dVqlgVdaC+uFj7fb/8APT9DVezv7cWifvPXsfWrG+q9m/8Aoifj/OmsErbE/XBIr+3+13H7z+72PpSXV/b/ALn95/y1XsadE/8ApVx/wH+VJct/qf8ArqtaLBpPYl4olN9b/wDPT9DVS1vYPsqfvPXsfWrm6q9of9FT8f511UqCj0OWpXUomj4VkSXV9UZDkeTbjp7y1oeKuNCP/X1bf+j46o+Gf+Qzqn/XG3/nLV3xZxoDf9fVt/6Pjr6iKvhLf3X+R69B2wyfkZ26pNM1JrTwldzwIfMt558mWNgi5mf5s8blAOTg9jyKq7q2vDH/ACAk/wCvi4/9HPXlZJR9nUm/IxwVXnk0Z48TkadY3L3unAvqBtJT0DjeVDJ8/wApK7XGd3DD1zVUeK75oNXZ2sLd7VyIVbBIAeRNr7pEAc7NwGRwe9dlRX0R6JBZTtc2NvcMNrSxK5GCMEjPQ4P581PRRQAVjatok2oXaXEF2kDAR7hJCZATHIJExhhj5hz6j061s0UAc5L4Yke485L5F8qV5rcGDO0vMszb/m+cblAGNuB69aTXdNjtPh5qFk0kji306Ub1Yx7isZ9D09unaukrE8YsV8D6+w6jTbg/+Q2ppXdiZfC7HJg2zHAAP/AjUdskH735P+Wh/iNef2muXRIO+P8AWtO01q5+f54/vnua9WtwzTi7RR+cSr4yne7uddJHAbqD5P738R9KW7ig+yv8np3PrXO/2vcNPCSYeN3r6VYn1Odrdx+57evrXDPIOXZErMK6audD5MH9z9TUNtDB+++T/loe5qguoTH/AJ5frTba8mPm8xf6w+tc8ska6AsznZ3Zo+RB9u+5/wAs/U+tF3BB9lf5PTufWqIuZze/eh/1fv6066nuDbP80Pb19ayeUtdC1mUrrU0/Ig/uf+PGoba3g/ffJ/y0Pc1F5tz6w/rUNtJc/vf9T/rD61m8st0Gsxdty39ng+3fc/5Zep9aLy3g+yv8np3PrVXzLn7b/wAsf9X7+tF3Jc/Zn/1Pb19an+zvItZg7rU0fs8H9z9TUFpbwfZU+T1/iPrTPMuv+mP61DaPc/Zk5h7+vrUvLvIax7tuWfs8H277n/LP1PrRc28H7n5P+Wg7mq++6+2/8sf9X7+tFw91+6/1P+sHrS/s8pY53Wpf+zwf88//AB41BZ28H2VPk9f4j60m+6/6Y/rUNo919lT/AFPf19an6gCxrtuT/Z4Pt33P+WX94+tLc20H7n5P+Wq9zVffdfbf+WP+r9/Wi4e6/df6n/WD1pfUSljHdal/7NB/c/U1BZ20H2VPk9e59aTfdf8ATH9ahtHuvsqf6nv6+tL6kCxjtuTxW0H2qf5P7vc+lFzbQfufk/5ar3NQRPdfap/9T/D6+lLcNdfuv9T/AKwetS8JYtYp33Lv2WD/AJ5/qags7WA2qfJ69z604Ndf9Mf1qOyN19lT/U9/X1rKdDlRcK7ktzd8KRpFq2qKgwPJtz195av+LuPDz/8AXzbf+j46oeFPM/tbVPM258m3+79ZaveMOPDj/wDXzbf+j469Cir0EvI+qpu2Cv8A3X+Rib63/C3/ACAU/wCu9x/6Oeua3V0HhOV20IgwSKEnn2klcSfvX6YP4c4rPC0fZts8vJKvPUmvI3qK5fVvEt/ZRhhpr2PyO6m+MbCdlGREnlyHDt2z6cA9rt9rq2us2VtFc2kkdwkw8kH960kfZTu9QVIx1HWu0+jNuiuJg8XX0mjQ3ckmnic3SRGOMbldWVGIUmQEld+DgMc8BT37agAoqOYStC4gdElI+RnQsoPqQCM/mK564s/EB1O3eW4E6h4islpm3ijAkBl8xDI2/cnA64PYdaAOlrB8bnHgHxGfTS7n/wBFNVa8ttQvtSumFrrFpbvZ7N0F4gZ3YL91TIVRlxjIHJJ5x96p4iTULf4T61HcIXnXTbpW8+bL+XsfBZhuBfbjOOCc4IFXT+NCex8322r4A+9Wla6wRu4f7x7Vx8chQ8Vbt7vBOT3r9dq4RSasfPVsFB3sjtU1g+ZHw/ftVxtYzCw2v+VcZHe/OnPrVz7dmM81wVMFozzamAjdaHaprH+zJ+VSW+r/AOs+ST757Vykd971PBfff5/iNcs8HqcM8DGz0OqXVv8ASs7JPuenvUs2rZgYbJO3b3rl0vv3+c/wf1qeS+zC3Nc8sHoznlgldaHWLq3+xJ+VJb6qR5nyS/6w9q55b73p0F79/n+M1jLBq+xzvBqz0OjXVT9rzsl/1f8Ad96fcaqTbsNkvb+H3rnlvf8ASc5/g/rUk15mBufT+dYywSs9CPqtpLQ6Yat/0yk/75qO01TFumYpO/8AD71ii896La7/AHC8+v8AOs3gI32Mvq7UdjcGpj7ZnyZf9X/d96dPqSnyv3Mv+sH8NYgu/wDSuv8AB/Wny3f+r5/jFZPARD2TTR0A1Jf+eU3/AHzUdpqC/Zk/dS9/4fessXfvTLW7/cLz6/zqHl8bkcklE2Bfr9sz5Uv+r/u+9LPfKfK/dS/6wfw1ki6/0vr/AAf1p811/quf+WgrJ5eg95NaG0L0f88pf++ajtLwfZk/dS9/4feqIuvem2l1/o6c+v8AOoll5PNJR2NGK7/0mf8Acy/w/wAPtS3F1/qv3Mv+sH8NUYbo/aZuf7v8qsS3OfKz/wA9BXNPAtItVmpLQvi6/wCmM3/fNNsbn/Rk/czd/wCH3p6PkZpbA/6Mn4/zryMXQ5Ys7cLWUjc8KP5mraodrL+5t+GGO8tXfGPHhuT/AK+bb/0fHVXwv/yF9T/642/85aseNTt8MSn0uLb/ANHx1nhVeEUfep2y9v8Auv8AI5vdXVeE/wDkXov+u0//AKOeuL88V2fhE58NwH1ln/8ARz13V6Hs0mfOcLVHOtU9F+Zt0UUVzH2oUUUUAFFFFABWB45/5J/4k/7Bd1/6Kat+sDxz/wAk/wDEn/YKuv8A0U1VD4kB8bjFOXHP1qIE+1KpPPTrX6zTrRTWhwNEwOGGKnEwC1Uy2R0pxJx2rp56clK6/qxDhc1I5l/yanhmTn6+tYwZh0IpySyLnkdacqUG1ZGEqCZvJMnm/wDAfWp2mTyzx+tc+ty+/t0qf7S+w8isJYePK9DCWGdzolmT/Jp8MsfzcfxHvWEt0/qtSRXMnzcr1rKeHjdaHNLDOxvrLH5/T+H196kkkj8o8frWCtzJ5vVfu1K9zJ5Z5WsJYdWehi8M7o6FZI/T9aW3ePyl4/WsRbmX1WnQ3MnljlKiWHV9jB4Z2N1Xj+0dP4PU+tPkeL5OP4x3rDW5l87qn3f61I1zL8nKfeHrWLoKz0M3h3fc6ANF/d/U0lu0Xkrx+p9ax1uZfVP1pYLmXyl5T9al0FfYxeHdtzaDRfaen8HqfWnymL938v8AGO5rFFzL5/VPu+/rUj3M3ycp94etZOgrPQh4d3WpvAxf3f1NFt5XkL8vr3PrWQtzN6p+tLbXE3kryn61MqCvsYug7bm3D5Xny/L6dz6VZYQ5i+X/AJaDuawobibz5OY+3r6VaM85MfMf3x61zzoq2xhOi77nVW6QHHy/qauafawG2T5PXufWsK2muOMGL9a2tMa5NpGR5P6+tfPZrRUYXSOfD80Z2udR4TjSLVdUVBgeTbnr7y1P48bb4QuD6T23/o+Oq/hLzP7V1TzNu7ybf7v1kqT4htt8FXZ9Jrb/ANHx15GBV3Bea/M/Roq+Xtf3X+RwX2mvQvBUryeGIh5Dptkm2s5GH/ePyMEkDtyBXk/2n3r1rwOd3g+yPq0p/wDIr17+cU+WEfU+e4Xpclap6L8yDVr3xJHGCtqludjmP7C7XRklA+RH3RDYhOct+q959UvtYOtadFp9vcLYs8iTzeUpGdkoG4MQwVWVDkD5twwezdDRXgH2ZxEGoeJvKj8z7Uf3jiJxbZMxAh2rJmJSiEmbnauAo57nt6KKAI5plgheVw5VBuIRC7fgqgk/QCsR/E6NeJFb2cxiDRLM9xHJbuvmyCNNqOgLc9TwAPU8Vv1DLZ2s88U81tDJNDnypHQFkz12k8j8KAMbUfENzpuozQT6Zi3MebWfzxid/l+UqASgBY5J7DODVDxXqX2r4YeIriWFo3Nje27JGGkAdRJGeQPu5UnJA464roV0jTEZmXTrQF4vJYiBfmjwBsPH3cADHTgVk+M4Ibb4ceIoYIkiiTSroKkahVX903QCqj8SA+Nt3saUN14NNzSg9a/SFN3Wpx2HbuRwaUtx0NMzyKUnit41ZWeorEgf2NKr9eD1pmaVT1+tdca0rrUmw/d83Q9KUv8AL3pmfm/ClJ4rdVW4vUViYTMPWnJcEZ4PWoc0A9frW102tSHFFtbr5+h6VIbr5Dwao5+b8KUsQvWk6d4vUh04mot37NT4rvCDhvyrMEpFOjmwoqZUXzLUydFWNVbv970bp6VIbv7vDfeHaspZvn69qkM33ee9ZSouzMnRXY11u/8AZb8qdDd4jHyt+VZiz+9Oin+Qc1EqL5jJ0VbY1Fu/333W+76VI139z5W+8O1Zgn/eZz2qRpvu89xWUqLszJ0VfY11vP8AZf8AKnW93+7Hyv8AlWas3vT4Jv3Y5rOVF3MZUVbY1obv96/yP27VcF3zH8j/AHx2rGhm/ev+FXVl5j5/iFctSk+VnLVpK+x1Fpd9P3cn/fNdJo9zmyj/AHUp6/w+9cnZy5xXWaE2bCP8f5185nFP90zxZpRmtOp1fhJ/M1XVDtZf3Nvwwx3lo+JTbfAd8fSW2/8AR8dP8Lf8hbU/+uNv/OWoPik234d6ifSS3/8AR8dfM5erzprzX5n6JRV8Cl/d/Q8d+017d4CO7wVp59fMP/kRq+evtPua9++HTzt4GsC0cYXa/lkSElvnbrxxz9a+q4ip8tKHr+h5mSUuSpP0Ororl9WvfEkcYK2qW52OY/sLtdGSUD5EfdENiE5y36r3fqF/qv8Awk1lFZrdCwO1ZcW52k73VuSh6YByWQYwR5mcV8mfRnS0Vytvcate6TaCK61KO5OoyQtJLZrGWhErEM4aMYHlLgEAZZgOtdVQAUUVTudUs7S7itZpSssuNoCMQMnau4gYXJ4GSMngUAXKwPHP/JPvEn/YKuv/AEU1XL7xBp2mTSxXks0Rjhe4Zvs8hXy0ALEMFwcZHAOecVmeLbyDUPhlr93avvgl0i5ZG2kZHlNzg81UfiQHxvgUADmk59aBn1r726utDlFwMilIGKbzkc0pzjrVqSs9BDsClAHNN59aBnnmt1NXWgh2BupTjFN5z1oOcda1U1yvQViTAoAHP1pvPrQM8810KorrT8hWH4G78KUgYpnO7r2pTnHWtVUXK9PyFYkwKFAx0pnPrQucda6VUXMvd/IVh+0bvwpT2+tM53de1Kc8c96vmi0/d/IViYMPSnRsu0VDz60i5x1q3GLkvd/IlxLasu/p2qQsvy/WqILBvvdvSnb24+bv6Vm6as/d/Ih0zTVl9KkhK7Bx+tZ6u3979Kkhdto+b9KiVJc3w/kYSp6GrDs8xuPTvV5AhKcfxDvWNC7eY3z+nar8bPuT95/EO1cdSmuV+7+RyVYeZ1VkkZx8v6muz8P28LWKZT9T61wlj5nH73/x2u38OJMbFcT4/wCAD1r5nOYL2L0Pm8WmpLXqjuPCUaRarqioMDybc9feSq/xZbb8NdUPo9v/AOj46seElddV1MO+8+Tb84x3kqn8YG2/C3Vz6Nb/APo+Ovk8sV61Jf3l+Z99hVfBxXkfPH2n3r6V+G53fD7SD6xsf/H2r5T+019VfDI5+HGiH1hP/obV9pxZT5aNP1f5GWApckmdZRRRXwx6YUUUUARzQRXMLwzxJLE42ujqGVh6EHrWUPC+lxXCS2kC2SbkeSG1jSNJSjb0LALnIYZ4Iz0ORWzRQBQOmAzm4N3cG5+zi3WY7MoM5LAbcbmOM8Y+VeBisLxdpkNn8MfENsskzhdPvJi5fazuyu7E7cDBZjxjHbGK6ysDx1/yT7xL/wBgq6/9FNVR3QHxdn2NAP1ooHevubu5yhnkcUpPHQ0ncUp6U7uz1AXPtQD7GigVqm7rUQZ+boaUnjoaTvSnpVpuz1AXd7GgHrwaKB3rdSd1qIXdz0NKW46Gk/i/Cg9K0UnyvUQ7d7GlVuOhpKF6V0KT5lqIXd83Q9KUt04PWkz81Ke31rRSdnqIdu9jQrcdDSUL0rpUpcy1EO3fN0PSgnpwetJn5vwpSen1rRNtPUQ4Pjsakil4HBqOhDgCtHFuS1JaTRdhl+duD27VoRTcp8rfeHasqBvnNaMLfMn+8K5KsHyvX+rnJWijp7CfGPkf8q7vwzcYsl/dSn6L71wVg3SvQfCx/wBCX/Pevls6i/YyPlcwsmtOqO38JP5mqamdrL+5txhhjvJVD4yf8kq1n62//o+OtHwr/wAhbU/+uNv/ADlrO+Mv/JKNa+tv/wCj46+Myx2rUn/eX5n3mC/3SHofK+6vrD4Wmdvhvo4ZI1UQny2DFifmPJGBjntk18l7q+uPhd/yTPQv+uB/9CavtOLK3tKNNeb/ACLoKzZe1C31xY7Z2kF5Il1CwFghttqBvn375iHBXjH6ek17/aDa1ZT2tndbFSZJN1wqxDshZQ/OSM5CkgHtyK26K+GOk4tNP8XLbxI7ozJqKTSEXzbpE3xk4O37mPNynHbGeh7SiigAoqOYStC4gdElI+RnQsoPqQCM/mKxJINYTVdLkmLXQWVxPJaDyIljKEDfG0p3Hdg5AP8AiAb9YHjr/kn3iX/sFXX/AKKaqV/a69c3d1NbJdw+YpZAblQojMKjygobAl8zcd4GB/e7VNeaRqWpeBdS0ku0VxdxXMK/aW81ljcuFUtu67Sozk496admB8aYFAAr6F/4Zy03/oYbv/vwv+NH/DOWm/8AQw3f/fhf8a+l/tTDef3GPIz57wMiggYr6E/4Zy03/oYbv/vwv+NRx/s6WTPKJNeulUPiMiFTuXaOTzxzkfhT/tXDa7/cHJI8AwKABzX0J/wznpv/AEMN3/34X/Gj/hnPTf8AoYbv/vwv+NX/AGthb9fuF7OR894GaUgYr38/s62QuEUa7cmIoxZ/JXIbIwMZ7gt+XvUn/DOmm/8AQw3f/gOv+NV/a+Es9/uD2cj59wKABzX0F/wzppv/AEMN3/4Dr/jUc/7O1klvK0Gu3MkoUlEaFQGbHAJzxzVrOcJfr9wvZyPAcDP4UpAxX0F/wzpp2f8AkYbv/wAB1/xo/wCGddO/6GG7/wDAdf8AGrWdYOzWv3B7OR8/YFAAxX0D/wAM7ad/0MN1/wCA6/41HB+zvZmM+brtyjb2AAhU/LuO09e4wfbNarPcFfr9wvZSPA8Dd+FLgcfWvoD/AIZ207Of+Eguv/Adf8aP+GdtO/6GC6/8B1/xq1n2Bt1+4PZSPAMChQMV74f2ebMXCKNduTEUYs/krkNkYGM9wW/L3p4/Z308D/kYLr/wHX/GtVxDgL31+4XspHgOBu/ClIHH1r37/hnjT85/4SC6/wDAdf8AGo5v2erRUBh124dt6ggwqPl3DcevYZPvirXEeX2e/wBwvYzPB8ChVG3pXvv/AAzzp/8A0MF1/wCA6/40D9nrTwP+Rguv/Adf8a3XE2W36/cL2MzweFRvPFaEKLuTj+IV7Uv7PlghJHiC5/8AAdf8adb/AAHg27pdcmRw7YAgU5UMdp69xg/jWU+JMucbK/3GVTDVJbHnFhEhx8v616B4XtoWs1yn6n1rZh+DUcH3ddkP1th/8VWxYfD+506IRw6yhA/v2mf/AGevDzLNcLiKbjTvf0PBxmTYurbktv3LvhKNItU1NUGB5NuevvLWf8Z/+ST619bf/wBHx1v6Fot3pWo3slxcx3Ec0UKoyRbMFTJkEbj/AHlp3i7w3F4u8L3mhzXD28dzszKigldrq/Q/7uK+Ywr9lKLl0f6n0+Gpyp0IwlukfF26vr34W/8AJMtB/wCvc/8AoTVwX/DOem/9DDd/9+F/xr1HwtoMnhnw9b6OLtLiK1TZDJ5JQ4yT83zHJ57Yr2s1zCnioRjDozSEHE26K57ULfXFjtnaQXkiXULAWCG22oG+ffvmIcFeMfp6O1CDWrrXtLmgQx6bFK4nT7SY2IKSqWYLkMv+rKjOQTyAeV8Q0N+iuWi0/UrRNOkWPUZXOoPLMn24uIoSGVVbfJ8wA2HA3c7iK6mgAooooAhuoJLi3aOK6mtXOMSwhCw/77Vh+lY6aVf2uuJex+Rd7bOSFri4kEc0rFgyhgkQXaMYz2z0PfeooA5GbwvePdXEuyykVpmkZHc4uwZvMCy/KcbF+VfvfgOK6PSrWSx0iytJXDyQQJGzDoSFAPWrdFAGBqehXd1tY3z38Sy+YbC9ESwOOeCVi3cZBGd3QZ9Q2Xw+95pel2F7Haz/AGOSHfPJ87sqKpJGRwWdQCM/dyc54roaKAMTQtFbSr7UpUggt4LqQSCKKZpcybnLuSwBG7cvyjIG3jrVjVNLu79JRb6xeWe+MoEiWIrnB5yULfkw9sVp0UAczF4el/sW80+fTrB7e5kEn2SO5eOOPAjAUMEzglGcnA5OOck1Hpnhq9s9R0+4ma2YWyIDIrHeqiJ08lfl5TLBs5GSOldVRQBDdQSXFu0cV1NaucYlhCFh/wB9qw/SsBNBv7fWlvhPHeOjGRbi4YJMw8pkEJ2IAE3EPkdx90nmulooA5K78MSX0uqtd6Vpskd1Inlxx3TwhwjFg8hWPO/JyTz6diW6Wxhlt9PtoJ5BLNHEqSSKu0MwABIHbJ7VYooAwLnRrtvE8OpxLbOiyKxd3KyBNjI0Ywp+XLb+vLDGO4xn8HXzWRh22LHaqMjO22ZxHIpuG+X/AFhLqSOfu/e6EdxRQA1FKRqpYsQACx7+9Z+t6fJqNlHEkcEwSZZGguDiOUD+Fjg/XoeQK0qKAOXj0K/36KLi20+UWEASS485hM7KCowdmQpB3EZ6nGSAd1/wzox0PTpLQQw28Hm7obeGVpFiTaoxvYAtyGPI43Y7Vs0UARXUTT2k0KSGJ5EZVkXqpIxkfSuTPhWeXQ7nT7jTtOeN5kmhhS7kjSJhGEJBEfqu7pyXOemT2NFAHPaVoV3Za699cSW8mYmVp0yJJiwixuXGAF2Njk/e7c56GiigDnrbRLmHxFNePHayW8wmSVy58yVX2FQy7cHZs2gZ6MTx0NJPCAiXSnis7JLq2uWuHmSRlEW+RXkREC4cFVCAtggAHrXXUUAFZOv6fcaha2y2sFrLNDdwzhrhymwJIrEqQrHJAI7deta1FAHJ6j4bv9QurmfNpFJOhJlEjFhmDy/IztGY93z545/hzzW3otjJp2neRIkUZMsjiKE5SMM5YKvA4APoK0aKACiiigAooooAKKKKAP/Z", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "scheme = solutionScheme(form==0, solver=\"gmres\",\n", " parameters={\"newton.linear.preconditioning.method\":\"jacobi\",\n", " # \"newton.verbose\": True,\n", " }\n", " )\n", "uh = space.interpolate(0, name=\"solution\")\n", "info = scheme.solve(target=uh)\n", "\n", "if parallel:\n", " gridView.writeVTK(\"u_h\", pointdata=[uh])\n", "else:\n", " uh.plot()" ] }, { "cell_type": "markdown", "id": "9a46fb5d", "metadata": {}, "source": [ "So far the example was not really advection dominated so we now\n", "repeat the experiment but set $\\varepsilon=1e-5$" ] }, { "cell_type": "code", "execution_count": 4, "id": "4bdda44b", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:32:55.063483Z", "iopub.status.busy": "2024-02-29T12:32:55.063250Z", "iopub.status.idle": "2024-02-29T12:32:55.534345Z", "shell.execute_reply": "2024-02-29T12:32:55.533278Z" } }, "outputs": [ { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrAS0DASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+ql/frp8cMjwyyJJPHCWj2/IXYIpOSOMsOmTz0q3VTUNNttTgSG6EpRJFlAjmeM7lOVOVIJwQDjpxQBRvfElpYXE8c0NwVhJRpVVSpkEfm+WOc7tnPTHbOasJrdiNNe/u5VsYY3eOU3cip5bKSCCc7e3r0psnh/TJmZpYHkLR+W2+Z2DfJs3HJ5fbxv+9jvWb4ss4LH4c+IoLdWVBpt23zOWJJjckkkkkkk9aBpXdif/AITfwl/0NGif+DCL/wCKo/4Tfwl/0NGif+DCL/4qvkSkXv8AWuX6z5Huf2Mv5/w/4J9ef8Jv4S/6GjRP/BhF/wDFUf8ACb+Ev+ho0T/wYRf/ABVfIZ+8KG+6aPrHkH9jrX3/AMP+CfXn/Cb+Ev8AoaNE/wDBhF/8VR/wm/hL/oaNE/8ABhF/8VXyJSL3+tH1nyD+xl/P+H/BPrz/AITfwl/0NGif+DCL/wCKo/4Tfwl/0NGif+DCL/4qvkP+P8KG+6aPrHkH9jq3x/h/wT68/wCE38Jf9DRon/gwi/8AiqP+E38Jf9DRon/gwi/+Kr5EpF7/AFo+s+Qf2Mv5/wAP+CfXn/Cb+Ev+ho0T/wAGEX/xVH/Cb+Ev+ho0T/wYRf8AxVfIf8f4UN900fWPIP7HVvj/AA/4J9ef8Jv4S/6GjRP/AAYRf/FUf8Jv4S/6GjRP/BhF/wDFV8iUi/dFH1nyD+xlf4/w/wCCfXn/AAm/hL/oaNE/8GEX/wAVR/wm/hL/AKGjRP8AwYRf/FV8h/x/hQ3b60fWPIP7HX8/4f8ABPrz/hN/CX/Q0aJ/4MIv/iqP+E38Jf8AQ0aJ/wCDCL/4qvkSkX7oo+s+Qf2Mr/H+H/BPrz/hN/CX/Q0aJ/4MIv8A4qj/AITfwl/0NGif+DCL/wCKr5D/AI/wobt9aPrHkH9jr+f8P+CfXn/Cb+Ev+ho0T/wYRf8AxVH/AAm/hL/oaNE/8GEX/wAVXyJSL90UfWfIP7GV/j/D/gn15/wm/hL/AKGjRP8AwYRf/FUf8Jv4S/6GjRP/AAYRf/FV8hj7xobt9aPrHkH9jq1+f8P+CfXn/Cb+Ev8AoaNE/wDBhF/8VR/wm/hL/oaNE/8ABhF/8VXyJSL90UfWfIP7GV/j/D/gn2ZputaVrKyNpep2d8sZAc2s6yhCemdpOKnvb60020e7vrqC1to8b5p5AiLk4GWPA5IH415H+z//AMgzXP8ArtF/6C1dX8Xxn4X6uP8Aat//AEfHW6neHMeXPD8uI9jfrY2v+E38J/8AQ0aL/wCDCL/4qtSy1Kx1KAzafeW95GP4reVZB7cg4r438qvqb4Zj/i3Gij/piw/8fas6Vb2jtY68fl31SKlzXuac3iS1tvD82rzQXCJE0qG3wplLxsysoAbBOUbvjAySBmnXOvC2u7m3OnXjmC2NyXQxYZQQMDL5BPOMgZ2t6Un/AAi+jtYTWM1o1xbTFmaO5meYBmLFmXex2k725GDzVt9IsXjuI/JKrcRLDJ5bsh2KCAoIIKgZPTHU1ueWWYJVngjmUELIoYA9cEZqSora3jtLaO3h3eXGoVd7lzge5JJ/GpaACiioLm8tbJEe6uYYFdgimVwoZj2Gep9qAJ6wfG//ACIHiP8A7Bdz/wCimrSn1fTbWd4LjUbSGZEMjxyTKrKgGSxBPAA5z6Vk+MJ4br4deIZreWOaJ9LuirxsGVv3TdCKGOO6Pkfb7mkVevJ6075v7v60i7ueO/rXmH27SuIV+Ycmhl+U8mlO7cOP1obdtPH60xNKzDb7mkVevJ6075v7v60i7ueO/rSG0ribfn6npQy/KeTS/Nv6dvWht208frTFZWYbfc0ir15PWnfN/d/WkXdzx39aQ2lcTb8/U9KGX5TyaX5t/Tt60Nu2nj9aYrKzDb7mkVflHJp3zf3f1pF3bRx+tIdlcTb8/U9KGXpyetL82/p29aG3ccd/WmKysG33NIq/KOTTvm/u/rSLu2jj9aQ7K4m35+p6UMvTk9aX5t/Tt60Nu447+tMVlYNvuaRV+UcmnfN/d/WkXdtHH60h2VxAvzHk0MvTk9aUbtx4/Wht3HHf1pisrBt9zQq/KOTS/N/d/WhC20fL+tIdlc92/Z/GNM1zr/r4v/QWrrvi2M/DPVR/t2//AKPjrkfgBn+zNcyMfvov/QWrsfiqu74camvrJb/+j467FpR+R81JXzFL+8vzR82eV7V9N/DYY+Hmjj/pk3/obV86/ZTX0P8ADeaJvAenRRyJJJCrK6KwJU72OD6H61x4CfNJn0XFNJQo035v8jraKyTr0Q0VtTNndeWkjxyRfJvQq5RifmwQCp6E8Us/iCzt9Qns3EpeGEyM6qCpOVGwc5L/ADpxj+Id69M+KNWiq2n3sepaba30SOkdzCkyq+NwDAEA4JGeexNWaACsnXdFbWIQiXIgPlywuWj35SRdrYGRhumDzj0Oa1qKAMu40mS8lM1zLbzP9ja2Eb2+YsvjzCVLcq21RtzwAeTmsTxLpj6f8LtdsxdEGOwu5C0UYVcFXbYFO7agztAByABg119YPjf/AJEHxH/2C7n/ANFNSexUfiR8jbl9aRWXnnvUuKao6/WvLufduLuMLLuXmhmXaeaeR8y0OPlNO4nF2Y3cvrSKy8896lxTVHX60rjcXcZuXf17UMy7TzT8fP8AhQ4+U07i5XZjdy+tIrLzz3qXFNUdfrSuNxdxm5d/XtQzLtPNPx8/4UOPlNO4uV2Y3cvrSKy7RzUuKag+UUrj5XcZuXf17UMy8c96fj5/woYdPrTuLldmN3L60isu0c1LimoPlFK4+V3Gbl39e1DMvHPen4+f8KGHT607i5XZjdy+tIrLtHNS4pqD5RSuPldxgZdzc0My8c96eB8zUMOn1p3Fyuw3cvrSIw2jmpcUiD5BSuPldz3L4AHOma3j/ntF/wCgtXafE0bvAGoD1ltv/R8dcZ8A/wDkG63/ANdov/QWrt/iON3ga9HrNbf+j467JO2Hb8mfLVNMxX+Jfmjw77OfSvdvAAx4H01faQf+RGrx77P7V7H4EGPBtgPQy/8Aoxq8nJ5805eh9BxRU56NP1f5F0eGtMGnz2O26NtOxaRTezEklixIO/IyzEnB5709/DujS3JuZdMtZJmTYzyRhi4yp+bP3jlVOTk8Vp0V758WV7GyttNsYLKzhWG2gQRxxr0VR0qxRRQAUUVheJrnU7a2jOm+eGKyEGGHzS0gX92jDBwpPU8YwORmgDdrB8b/APIg+I/+wXc/+imqlqN7rK6jPJp893NbyWpZLc2BT7OdoO8Oy/O3X5OuT04IqHxDNfy/C7X3lVpJTZXq7rgeU5iAkCsVC/e24OMDPtSexUfiR8u7PdqRU68t1qTD+1NUPz0615Vz9AcFdaDSnzLy1Dp8p5anEPvXpQ4fYelFxOCs9BNnu1IqdeW61Jh/amqH56daLjcFdaDdnz9W6UOnynlqdh/M7dKHD7D0ouLkVnoJs92pFTry3WpMP7U1Q/PTrRcbgrrQbs+fq3Sh0+U8tTsP5nbpQ4fYelFxcis9BNnu1IifKOWqTD+1NQPsHSi4+RX2G7Pn6t0oZOnLdadh/M7dKGD8dOtFxcis9BNnu1IifKOWqTD+1NQPsHSi4+RX2G7Pn6t0oZOnLdadh/M7dKGD8dOtFxcis9BNnu1IifKOWqTD+1NQPsHSi4+RX2GhPmblqGTpy3WnAPvbpQwfjp1ouLkVthNnu1IifKOT+dSYf2pED7B92lcrkV9j3D4BjGna3/12i/8AQWru/iAN3g26HrPbf+j464X4CZ/s7W84/wBdF0/3WrvvHI3eE5x63Ft/6PjrsqP/AGWT/uv8j47EaZg/8SPMPs9eh+GNThsvAyyxFJ5bd5EaJZAMMZWChjzt6jnHTmuR+zmu88Jwxy+EoYZo1kjd51dHGQw81wQQeor5/h6fNVn6Hp55U56cPUhvPFj2dlbzGwDu8k6SosrERiJ9jtkIflz/ABMFUcbiuRVm5126t7HWJxZQM+nTmPYbggSr5aSAg7DhjvA24PPer50XSjGsZ0yy2I+9V8hcK2AMgY64AGfYU9NL0+OSeRLC1V55FlmZYVBkdTlWY45IIBBPSvqz5osRGRoUaVAkhUFlVtwU9wDxn60+iigAooooAKwvG3/IheIv+wXc/wDopq3awvGv/Ih+Iv8AsGXP/opqT2Kh8SPlDcvrTVK/Nz3qztpqL97614vMj9OdKV0QErvXn1ocrsPNTlfnX8aJF+Q0+ZEulKzIty+tNUr83PerO2mov3vrS5kU6Urogyvmde1Dldh5qfb+8/CldPkNPmRLpSsyHcvrTVK/Nz3q1spqJ9760uZFOlK6K+V8zr2ocrsPNWNn7z8KHT5DT5kS6UrMh3L601Cuwc1a2U1E+QUuZFeylcr5XzOvahivy896sbP3n4UOn3frT5kT7KVmQ7l9aahXYOatbKaifIKXMivZSuV8r5nXtQxX5ee9WNn7z8KHT7v1p8yJ9lKzIdy+tNQrsHNWtlNRPkFLmRXspXK4K7259KGK/Lz3qwE+dvwodPu/WnzK5PspcpDuX1pEZdg5qzspsafIKXMrFeylc9n+AxB0/W8f89Yv5NXoHjMbvDUg9bm2/wDR8dcF8Cxix1v/AK6xfyavQfFg3aAw9bq2/wDR8ddtV/7FJ/3X+R8Fj/dx8v8AEcn9n9q6zwmMeHYR6TT/APo56yPs9XNG1JbHwtPKsTu8FxMu1lZFZjOwHzEYxkjJGcc18twrPmrVPRfmXmNTnjE6aiubvPF8NhY2ctxFEk9xeG1MbThVXbN5TuGIGQCQQMZOR05IZ/wlVyb6Wx/sry7sSBY1mlZFZSJmyWKccQN93cORz1x9ueSdPRUFndJe2NvdxqypPGsihhyAwzz+dT0AFZOu602jwq6WwnPlyzOGk2YSNdzYODlumBxn1GK1qhubO1vERbq2hnVGDqJUDBW9RnoaAM2+1W+t7+a0tNOivJFtmmjRLnaxI+6r5XCBjkA5PQ8cHGV4k1WO8+GOuXVwEtXewu4WjaQHbKFdCmeMncCOOvbNbw0TSQ7ONLsgzxiJm+zpkoAAFPHQAAY9BWf4rt4bb4f67BbxRxRJplyEjjUKqjy26AdKT2Kh8SPmD7MfVvzpqWx+blvvHvW39im/up+tMis5vn+VPvH1r5r26tufrrjTutDHa2O9OW796JLYhDy351sNZzebH8qd/Wiazm8pvlT9aPbrTUlxp2lp/VjJ+zH1b86alsfm5b7x71ufY5v7ifrUcdnN8/yp94+tHt/Mpxp3Whj/AGY+b1b7vrRJbEIeW/Otj7HN5/3U+77+tE1nN5TfKn60e31WpLjT5XoZP2Y+rfnTUtj83LfePetz7HN/cT9ajjs5vn+VPvH1o9v5lONO60Mf7MfN6t931oktiEPLfnWx9jm8/wC6n3ff1oms5vKb5U/Wj2+q1JcafK9DJ+zH1b86bHbEoOW/Otz7HN/cT9ajhs5vKX5U/Wj26tuVy0+ZaGP9mPm9W+760PbH5eW+8O9bH2Obz/up9339aJLOb5PlT7w9aPb+ZPLTs9DJ+zH1b86bHbEoOW/Otz7HN/cT9ajhs5vKX5U/Wj26tuVy0+ZaGP8AZj5vVvu+tD2x+XlvvDvWx9jm8/7qfd9/WiSzm+T5U+8PWj2/mTy07PQyfsx9W/Omx2xKDlvzrc+xzf3E/Wo4bObyl+VP1o9urblctPmWhjrbHe/Ldu9D2x+XlvvDvWwtnN5snyp29aJLOb5PlT7w9aPb67k8tPl2/q5k/Zj6t+dJHbExjl62/sc39xP1pkNpP5S/Kn60vb6blWp82x6d8DkMdnrQOf8AWRHn6NXoPif/AJAv/b1a/wDo+OuH+DkbxRayrgA7oTx9Hrt/FHGif9vVr/6Pjr15vmwMvOL/ACZ+ZZu0sdVa7lTatW9Aijm0F4ZY1kiee5VkcZDAzPkEdxVHdWh4a/5Aw/6+bj/0c9fMcK0uStUfkvzPPlW9poXotOsoIDDDZ28cRdXKJEoXcMYOAOo2rg+w9KY2kaY8c0b6daMk7+ZKpgUiRv7zDHJ9zVyivtyQAAGAMAUUUUAFFFQTXlrbzQwzXMMcsxIiR3AaQjrtB69R0oAnrG8X8+CteH/UOuP/AEW1Xn1TT42nV7+1VrfHnBplHlZ6bueM+9U/Ejxz+DtXeN1eN9PmKspyGBjPIPcVM/hY4uzTPGvs0Xqf++TUUNtF+85P3z/Ca9O/s+z9KgttPtD53H/LVq/NFi52fuM+q/1hpX+I84e2i8+Lk9/4T6UXFtEIG5Pb+E+tejy6fafa7cY/vfyovNPtBaOcen86axc7r3GH+sNLX3jz37NF6n/vk1FDbRfvOT98/wAJr0/+zrT0qvbafaHzuP8Alq1JYudn7jH/AKw0r/EecfZovtXU/c/un1ouLaIQNye38J9a9I/s+0+34x/yy/rSXmn2gtHOPT+dNYud17jF/rDSs/ePPfs0Xqf++TUUNtF+85P3z/Ca9P8A7OtPSq9tp9ofO4/5atSWLnZ+4x/6w0r/ABHnH2aL7V1P3P7p9aLi2iEDcnt/CfWvR/7PtPt+Mf8ALL+tF5p9oLRzj0/nTWLnde4xf6wUrP3jz37NF6n/AL5NRW9tEYF5Pf8AhPrXp/8AZ1p6VXs9PtDaIcev86X1udvgY/8AWGle/MecfZovtXU/c/un1omtov3fJ++P4TXo/wDZ9p9vxj/ll/Wi50+0Hk8f8tVp/W53XuMX+sFK3xHnv2aL1P8A3yait7aIwLye/wDCfWvTv7OtPSoLPT7Q2iHHr/Ol9bnb4GP/AFhpXvzHnH2aL7V1P3P7p9aJraL93yfvj+E16P8A2fafb8Y/5Zf1oudPtB5PH/LVaf1ud17jF/rBSt8R579mi9T/AN8more2iMC8nv8Awn1r07+zrT0qCz0+0Nohx6/zpfW52+Bj/wBYaV78x5wltF58vJ7fwn0omtov3fJ++P4TXo8Wn2n2u4GP7v8AKi50+0Hk8f8ALVaf1ufN8DF/rBStbmPPfs0Xqf8Avk1Hb28XkLz6/wAJ9a9N/s609Kgs9PtDaJx6/wA6X1udvgY/9YaV/iIPhkipNq4U8Yh7f79dP4q40I/9fVt/6Pjqn4ahjg1rVFi+75Nuf1lq34s40Bv+vq2/9Hx19zhG55cm1vFnzGOrqtUnVXW7M7dUmmak1p4Su54EPmW88+TLGwRczP8ANnjcoBycHseRVXdW14Y/5ASf9fFx/wCjnrhySj7OpN+R5OCq88mjPHicjTrG5e904F9QNpKegcbyoZPn+UldrjO7hh65qqPFd80Grs7WFu9q5EKtgkAPIm190iAOdm4DI4Peuyor6I9Egsp2ubG3uGG1pYlcjBGCRnocH8+anoooAKxtW0SbULtLiC7SBgI9wkhMgJjkEiYwwx8w59R6da2aKAOcl8MSPcecl8i+VK81uDBnaXmWZt/zfONygDG3A9etJrumx2nw81CyaSRxb6dIN6sY9xWM+h6e3TtXSVieMWKeB9fYdRptwf8AyG1NK7sTL4XY5Pda+g/76NQW7W373gf6w/xGvNP+Eju/7yfrTIvEV2N/zJyxPeve/wBUqXY/P0sfb4j06Vrb7VBwP4v4j6UXbW32V+B2/iPrXmL+IrsyxncnGfWibxFdtEw3J+tL/VOl22Glj7r3j1Xda+g/76NQW7W373gf6w/xGvNP+Eju/wC8n60yLxFdjf8AMnLE96f+qVLsJLH2+I9O3W327oP9X/ePrRdtbfZX4Hb+I+teY/8ACRXf2jO5Pu47+tE3iK7aJhuT9aX+qdK17Dtj7r3j1Xda+g/76NQW7W373gf6w/xGvNP+Eju/7yfrTIvEV2N/zJyxPen/AKpUuwksfb4j07dbfbug/wBX/ePrRdtbfZX4Hb+I+teY/wDCRXf2jO5Pu47+tE3iK7aJhuT9aX+qdK17Dtj7r3j1Xdbeg/76NQWjW32VOB3/AIj615p/wkd3/eT9aZD4iu1iUbk/Wn/qnS2sTbH2+I9O3W327oP9X/ePrRcNbfuuB/rB/Ea8x/4SK7+0Z3J93Hf1ol8RXZ2fMnDA96X+qdLsVbH3XvHqu629B/30agtGtvsqcDv/ABH1rzT/AISO7/vJ+tMh8RXaxKNyfrT/ANU6W1ibY+3xHp262+3dB/q/7x9aLhrb91wP9YP4jXmP/CRXf2jO5Pu47+tEviK7Oz5k4YHvS/1TpdirY+6949V3W3oP++jUFo9r9lT7vf8AiPrXmn/CR3f95P1pkPiK7WFRuT9af+qdLaxNsfb4j06J7X7VP93+H+I+lFw9r+6+7/rB/Ea8xTxFdCWQ7k5x60S+Iro7PmThge9L/VOla9irY+/xHqm+19F/76P+NQ2jW32VOB3/AIj615r/AMJHd/3k/Wo4fEN2sSjen5mn/qnS2sTbH2+I9u8HmM6pqvl9PKg757yVpeLuPDz/APXzbf8Ao+OuP+Ed7JfLrMspBYGFePo/+Ndf4w48OP8A9fNt/wCj468TFYZYeUqC6aH2OGc1gE578rMTfW/4W/5AKf8AXe4/9HPXNbq6DwnK7aEQYJFCTz7SSuJP3r9MH8OcV52Fo+zbZ52SVeepNeRvUVy+reJb+yjDDTXsfkd1N8Y2E7KMiJPLkOHbtn04B7Xb7XVtdZsraK5tJI7hJh5IP71pI+ynd6gqRjqOtdp9GbdFcTB4uvpNGhu5JNPE5ukiMcY3K6sqMQpMgJK78HAY54Cnv21ABRUcwlaFxA6JKR8jOhZQfUgEZ/MVz1xZ+IDqdu8twJ1DxFZLTNvFGBIDL5iGRt+5OB1wew60AdLWD43OPAPiM+ml3P8A6Kaq15bahfaldMLXWLS3ez2boLxAzuwX7qmQqjLjGQOSTzj71TxEmoW/wn1uO4QvOum3St582X8vY+CzDcC+3GccE5wQKun8aE9j5f8A7TP+1+VNXUiM/e6+lZ+6kDdfrX7G4wutTyfYQ7GidSO9T83ftQ+pEoR835VnFuRQzcVLjC0tf6sP2EOxp/2mf9r8qaupEZ+919Kz91IG6/WqcYXWovYQ7Gj/AGkfMz83T0ofUiUI+b8qzt3zfhQzcVLjDleo/YQ7Gn/aZ/2vypq6kRn73X0rP3Ugbr9apxhdai9hDsaP9pHzM/N09KH1IlCPm/Ks7d834UM3FS4w5XqP2EOxp/2mf9r8qampEIB835Vn7qRW4quWHMtRewh2NH+0j5mfm6elDakTj73X0rO3fN+FBbp9alxhZ6j9hDsaf9pn/a/KmpqRCAfN+VZ+6kVuKrlhzLUXsIdjR/tI+Zn5unpQ2pE4+919Kzt3zfhQW6fWpcYWeo/YQ7Gn/aZ/2vypqakQgHzflWfupFbiq5Ycy1F7CHY0RqR3sfm7dqG1InH3uvpWcG5NBbp9anlhy7/1cfsIdjT/ALTP+1+VImpkIB835VnbqFPFVyQ5lqL2EOx9B/AOf7Rp2tvzxNEOf91q9D8Y8eG5P+vm2/8AR8debfs8nOk65/13i/8AQTXo/jU7fDEp9Li2/wDR8dflubxvmNRL+Y6qyUcHO38r/I5vdXVeE/8AkXov+u0//o564vzxXZ+ETnw3AfWWf/0c9cleh7NJny/C1RzrVPRfmbdFFFcx9qFFFFABRRRQAVgeOf8Akn/iT/sF3X/opq36wPHP/JP/ABJ/2Crr/wBFNVQ+JAfG3FIMc/Wm5PtQCeelfqDqK60OKw44yKDjFNJOR0oJOO1J1FZ6BYfxSDHP1puT7UAnnpTdRXWgWHcbvwoOMU3J3dqCTjtSdRWegWH8Ugxz9abk+1AJ56U3UV1oFh3G78KDjFNyd3agk47UnUVnoFh/FIMYpuT7UAnHan7RX2Cw7jd+FBxx9abk7u1BJ46UvaKz0Cw/ikGMU3J9qATjtT9or7BYdxu/Cg44+tNyd3agk8dKXtFZ6BYfxSDGKbk+1AJx2p+0V9gsOGMmg44+tNBOT0oJPHSl7RcuwWH8Ugxim5PtQCcdqftFfYLH0F+zv/yCNc/6+Iv/AEE16L48bb4QuD6T23/o+OvOf2df+QPrn/XxF/6Ca9C+IbbfBV2fSa2/9Hx18DjfezGX+IrEK+Fmv7r/ACOC+016F4KleTwxEPIdNsk21nIw/wC8fkYJIHbkCvJ/tPvXrXgc7vB9kfVpT/5FeuvOKfLCPqfL8L0uStU9F+ZBq174kjjBW1S3OxzH9hdroySgfIj7ohsQnOW/Ve8+qX2sHWtOi0+3uFsWeRJ5vKUjOyUDcGIYKrKhyB824YPZuhorwD7M4iDUPE3lR+Z9qP7xxE4tsmYgQ7VkzEpRCTNztXAUc9z29FFAEc0ywQvK4cqg3EIhdvwVQSfoBWI/idGvEit7OYxBolme4jkt3XzZBGm1HQFuep4AHqeK36hls7WeeKea2hkmhz5UjoCyZ67SeR+FAGNqPiG503UZoJ9MxbmPNrP54xO/y/KVAJQAsck9hnBqh4r1L7V8MPEVxLC0bmxvbdkjDSAOokjPIH3cqTkgcdcV0K6RpiMzLp1oC8XksRAvzR4A2Hj7uABjpwKyfGcENt8OPEUMESRRJpV0FSNQqj903QCqj8SA+Nt3saQHrwaKQd6/Q3J3WpyCk8jg0E8dDSdxQelJydnqA7d7GkB68GikHem5O61AXPPQ0E8dDSd6D0pOTs9QHbvY0gPXg0Ug71Tk7rUBc89DQTx0NJ3oPSpcnZ6gO3expAeOhopB0quZ33AXPPQ0E9ODSd6D2qeZ2eoDt3saQHjoaKQdKrmd9wFzz0NBPTg0neg9qnmdnqA7d7GkB46GikHSq5nfcBc8ng0E9ODSdzQe1TzO24Dt3saQNx0NFA6VXM77gfQX7Ohzo+u/9fEX/oJrv/iU23wHfH0ltv8A0fHXAfs6f8gbXf8Ar4i/9BNd38Um2/DvUT6SW/8A6Pjr42vrmP8A28vzRdVXoSXk/wAjx37TXt3gI7vBWnn18w/+RGr56+0+5r374dPO3gawLRxhdr+WRISW+duvHHP1r2OIqfLSh6/oeJklLkqT9Dq6K5fVr3xJHGCtqludjmP7C7XRklA+RH3RDYhOct+q936hf6r/AMJNZRWa3QsDtWXFudpO91bkoemAclkGMEeZnFfJn0Z0tFcrb3GrXuk2giutSjuTqMkLSS2axloRKxDOGjGB5S4BAGWYDrXVUAFFFU7nVLO0u4rWaUrLLjaAjEDJ2ruIGFyeBkjJ4FAFysDxz/yT7xJ/2Crr/wBFNVy+8Qadpk0sV5LNEY4XuGb7PIV8tACxDBcHGRwDnnFZni28g1D4Za/d2r74JdIuWRtpGR5Tc4PNVH4kB8bYFAA5o59aQZ9a+7bV9jlFwMiggYpOcjmg5x1pNqz0AXAoAHNHPrSDPrTbV9gFwM0EDFJznrQc460rqz0AXAoAHNHPrSDPrTbV9gFwM0EDFJznrQc460rqz0AXAoAGKOfWkGcdad1fYBcDNBA4pOc9aDn1pXVnoAuBQAMUc+tIM4607q+wC4GaCBxSc560HPrSurPQBcCgAYo59aQZx1p3V9gFwMmggcUnOTzQc+tK6tsAuBQAMUc+tIM4607q+wH0L+zn/wAgbXf+viP/ANBNdv8AFltvw11Q+j2//o+OuH/Zz/5A2u/9fEf/AKCa7T4wNt+Furn0a3/9Hx18vLXMl/iX5o1kr02vI+ePtPvX0r8Nzu+H2kH1jY/+PtXyn9pr6q+GRz8ONEPrCf8A0Nq+k4sp8tGn6v8AI4cBS5JM6yiiivhj0wooooAjmgiuYXhniSWJxtdHUMrD0IPWsoeF9LiuEltIFsk3I8kNrGkaSlG3oWAXOQwzwRnocitmigCgdMBnNwbu4Nz9nFusx2ZQZyWA243McZ4x8q8DFYXi7TIbP4Y+IbZZJnC6feTFy+1ndld2J24GCzHjGO2MV1lYHjr/AJJ94l/7BV1/6KaqjugPi7PsaAevBopBX3Dvc5RSeRwaCeOhpO4oPSpd7MYufY0A9eDRSCqd7iFzz0NBPHQ0neg9KnWzGLn2NAPXg0Ugqne4hc89DQTx0NJ3oPSp1sxi59jQDx0NFIOlVrcBc89DQT04NJ3oNTrYBc+xoB46GikHSq1uAueehoJ6cGk70Gp1sAufY0A8dDRSDpVa3AXPJ4NBPTg0nc0Gp1sAufY0A8dDRQOlVrcR9Cfs5/8AIG13/r4j/wDQTXY/GT/klWs/W3/9Hx1x37Of/IG13/r4j/8AQTXY/GX/AJJRrX1t/wD0fHXzMnbMk/7y/NG/2D5X3V9YfC0zt8N9HDJGqiE+WwYsT8x5IwMc9smvkvdX1x8Lv+SZ6F/1wP8A6E1fS8WVvaUaa83+RjQVmy9qFvrix2ztILyRLqFgLBDbbUDfPv3zEOCvGP09Jr3+0G1qyntbO62KkySbrhViHZCyh+ckZyFJAPbkVt0V8MdJxaaf4uW3iR3RmTUUmkIvm3SJvjJwdv3MeblOO2M9D2lFFABRUcwlaFxA6JKR8jOhZQfUgEZ/MViSQawmq6XJMWugsrieS0HkRLGUIG+NpTuO7ByAf8QDfrA8df8AJPvEv/YKuv8A0U1Ur+1165u7qa2S7h8xSyA3KhRGYVHlBQ2BL5m47wMD+92qa80jUtS8C6lpJdori7iuYV+0t5rLG5cKpbd12lRnJx7007MD40wKABzX0L/wzlpv/Qw3f/fhf8aP+GctN/6GG7/78L/jX0n9p4a/X7jHkZ89YGRQQMV9C/8ADOWm/wDQw3f/AH4X/Go4/wBnSxZ5RJr10qh8RkQqdy7RyeeOcj8KX9p4bX/IORnz9gUADmvoX/hnLTf+hhu/+/C/40f8M5ab/wBDDd/9+F/xp/2nhr9fuDkZ89YGaCBivoE/s6WIuEUa9dGIoxZ/JXIbIwMZ75b8vepP+GctN/6GG7/78L/jS/tPDWf+QcjPnrAoAHNfQv8Awzlpv/Qw3f8A34X/ABqOf9nSxS3laDXrqSUKSiNCoDNjgE545p/2nhr9fuDkZ8/YGaCBivoX/hnLTf8AoYbv/vwv+NH/AAzlpv8A0MN3/wB+F/xpf2nhrP8AyDkZ89YFAAxX0L/wzlpv/Qw3f/fhf8ajh/Z0sWjJl166Rt7AAQqfl3HaevcYP40/7Tw1+v3ByM+fsDNBA4r6F/4Zy03/AKGG7/78L/jR/wAM5ab/ANDDd/8Afhf8aX9p4a3/AAA5GfPWBQAMV9An9nSyFwijXroxFGLP5K5DZGBjPcFvy96k/wCGctN/6GG7/wC/C/40/wC08Nfr9wcjPnrAzQQOK+hf+GctN/6GG7/78L/jUc37OliqAxa9dO29QQYVHy7huPXsMn3xS/tPDW/4AcjPn7AoAGK+hf8AhnLTf+hhu/8Avwv+NH/DOWm/9DDd/wDfhf8AGn/aeGv1+4ORnz1gZNBA4r6F/wCGctN/6GG7/wC/C/41HD+zpYshM2vXSNvYACFT8oY7T17jB/Gl/aeGt/wA5GfP2BQAMV9C/wDDOWm/9DDd/wDfhf8AGj/hnLTf+hhu/wDvwv8AjT/tTDX/AOAHIw/Zy/5Auu/9fEf/AKCa7D4z/wDJJ9a+tv8A+j46m+H3gFfAI1G1hvWu7a68uRXdArBxuDDAPTG38zW34u8NxeLvC95oc1w9vHc7MyooJXa6v0P+7ivFnWi8V7VbXuaJe7Y+Lt1fXvwt/wCSZaD/ANe5/wDQmrgv+Gc9N/6GG7/78L/jXqPhbQZPDPh630cXaXEVqmyGTyShxkn5vmOTz2xXoZrmFPFQjGHRkwg4m3RXPahb64sds7SC8kS6hYCwQ221A3z798xDgrxj9PR2oQa1da9pc0CGPTYpXE6faTGxBSVSzBchl/1ZUZyCeQDyviGhv0Vy0Wn6laJp0ix6jK51B5Zk+3FxFCQyqrb5PmAGw4G7ncRXU0AFFFFAEN1BJcW7RxXU1q5xiWEIWH/fasP0rHTSr+11xL2PyLvbZyQtcXEgjmlYsGUMEiC7RjGe2eh771FAHIzeF7x7q4l2WUitM0jI7nF2DN5gWX5TjYvyr978BxXR6VayWOkWVpK4eSCBI2YdCQoB61booAwNT0K7utrG+e/iWXzDYXoiWBxzwSsW7jIIzu6DPqGy+H3vNL0uwvY7Wf7HJDvnk+d2VFUkjI4LOoBGfu5Oc8V0NFAGJoWitpV9qUqQQW8F1IJBFFM0uZNzl3JYAjduX5RkDbx1qxqml3d+kot9YvLPfGUCRLEVzg85KFvyYe2K06KAOZi8PS/2LeafPp1g9vcyCT7JHcvHHHgRgKGCZwSjOTgcnHOSaj0zw1e2eo6fcTNbMLZEBkVjvVRE6eSvy8plg2cjJHSuqooAhuoJLi3aOK6mtXOMSwhCw/77Vh+lYCaDf2+tLfCeO8dGMi3FwwSZh5TIITsQAJuIfI7j7pPNdLRQByV34YkvpdVa70rTZI7qRPLjjunhDhGLB5Csed+Tknn07Et0tjDLb6fbQTyCWaOJUkkVdoZgACQO2T2qxRQBgXOjXbeJ4dTiW2dFkVi7uVkCbGRoxhT8uW39eWGMdxjP4OvmsjDtsWO1UZGdtsziORTcN8v+sJdSRz9373QjuKKAGopSNVLFiAAWPf3rP1vT5NRso4kjgmCTLI0FwcRygfwscH69DyBWlRQBy8ehX+/RRcW2nyiwgCSXHnMJnZQVGDsyFIO4jPU4yQDuv+GdGOh6dJaCGG3g83dDbwytIsSbVGN7AFuQx5HG7HatmigCK6iae0mhSQxPIjKsi9VJGMj6VyZ8Kzy6Hc6fcadpzxvMk0MKXckaRMIwhIIj9V3dOS5z0yexooA57StCu7LXXvriS3kzEytOmRJMWEWNy4wAuxscn73bnPQ0UUAc9baJcw+Iprx47WS3mEySuXPmSq+wqGXbg7Nm0DPRieOhpJ4QES6U8VnZJdW1y1w8ySMoi3yK8iIgXDgqoQFsEAA9a66igArJ1/T7jULW2W1gtZZobuGcNcOU2BJFYlSFY5IBHbr1rWooA5PUfDd/qF1cz5tIpJ0JMokYsMweX5GdozHu+fPHP8Oea29FsZNO07yJEijJlkcRQnKRhnLBV4HAB9BWjRQAUUUUAFFFFABRRRQB/9k=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "eps.value = 1e-5 # could also use scheme.model.eps = 1e-5\n", "uh.interpolate(0)\n", "scheme.solve(target=uh)\n", "\n", "if parallel:\n", " gridView.writeVTK(\"u_h\", pointdata=[uh])\n", "else:\n", " uh.plot()\n", "\n", "# for parallel runs skip 1d part\n", "if parallel:\n", " sys.exit(0)" ] }, { "cell_type": "markdown", "id": "76b18afa", "metadata": { "lines_to_next_cell": 0 }, "source": [ "## A linear transport problem\n", ".. index:: Equations; Transport\n", "\n", ".. index:: Spaces; DG Legendre\n", "\n", ".. index:: Methods; Method of Lines\n", "\n", "We now use a DG space to approximate the solution to the transport problem\n", "$$ \\partial_t u + b\\cdot \\nabla u = f $$\n", "\n", "with a constant transport velocity vector $b$. We will solve the problem\n", "in one space dimension taking $b=\\frac{1}{2}$.\n", "\n", "\n", "We first setup the 1D grid and the Discontinuous Galerkin space." ] }, { "cell_type": "code", "execution_count": 5, "id": "829d4682", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:32:55.538972Z", "iopub.status.busy": "2024-02-29T12:32:55.538764Z", "iopub.status.idle": "2024-02-29T12:35:01.249388Z", "shell.execute_reply": "2024-02-29T12:35:01.248210Z" } }, "outputs": [], "source": [ "# remove later, shows same error: from dune.fem.space import dgonb as dgSpace\n", "N = 200\n", "gridView = leafGridView([0], [1], [N], overlap=1)\n", "space = dgSpace(gridView, order=1)" ] }, { "cell_type": "markdown", "id": "e1cfe2ee", "metadata": {}, "source": [ "The space discretization uses the same upwind flux discussed already for the\n", "Advection-Diffusion problem in a previous example. This can be written as\n", "an ODE\n", "$$ \\int_\\Omega \\partial_t u(t)\\varphi_i = F_i(u(t)) $$\n", "\n", "where $(\\varphi_i)_i$ is a basis of a DG space, i.e.,\n", "$\\varphi_i=\\varphi^E_l$ only has support on one element $E$ of the grid.\n", "This leads to the matrix vector problem\n", "$$ M\\frac{d}{dt} \\vec{u}(t) = \\vec{F}(\\vec{u}(t)) $$\n", "\n", "where $M$ is the mass matrix which in the case of a DG space is block\n", "diagonal - in fact with the `dune.fem.space.dglegendre` space the mass\n", "matrix is even diagonal. The ODE is thus of the form\n", "$$ \\frac{d}{dt} \\vec{u}(t) = M^{-1}\\vec{F}(\\vec{u}(t)) $$\n", "\n", "We discretize this ODE using an explicit two-step Runge-Kutta method with\n", "a fixed time-step $\\tau$. Given $\\vec{u}^n$ we define the solution at the\n", "next time step by\n", "$$ \\vec{u}^{n+1} := \\frac{1}{2}\\Big(\\vec{u}_1^{n+1} + \\vec{u}_2^{n+1} \\Big) $$\n", "\n", "with\n", "\\begin{align*}\n", "\\vec{u}_1^{n+1} &= \\vec{u}^n + \\tau M^{-1}\\vec{F}(\\vec{u}^n) \\\\\n", "\\vec{u}_2^{n+1} &= \\vec{u}_1^n + \\tau M^{-1}\\vec{F}(\\vec{u}_1^n)\n", "\\end{align*}\n", "This is __Heun's method__.\n", "\n", "We first setup the form for the right hand side $\\vec{F}$ which is very\n", "similar to the stationary case described above:" ] }, { "cell_type": "code", "execution_count": 6, "id": "fecee87d", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:35:01.254508Z", "iopub.status.busy": "2024-02-29T12:35:01.254269Z", "iopub.status.idle": "2024-02-29T12:35:01.264368Z", "shell.execute_reply": "2024-02-29T12:35:01.263118Z" } }, "outputs": [], "source": [ "u = TrialFunction(space)\n", "v = TestFunction(space)\n", "n = FacetNormal(space)\n", "x = SpatialCoordinate(space)\n", "\n", "# transport direction and upwind flux\n", "speed = 0.5\n", "b = as_vector([speed])\n", "hatb = (dot(b, n) + abs(dot(b, n)))/2.0\n", "# boundary values (for left boundary)\n", "dD = conditional(x[0]<0.01,1,0)\n", "g = 1\n", "\n", "# forms\n", "advInternal = dot(-b*u, grad(v)) * dx\n", "advSkeleton = jump(hatb*u)*jump(v) * dS\n", "advBnd = ( hatb*u + (dot(b,n)-hatb)*g )*v * dD*ds\n", "form = (advInternal + advSkeleton + advBnd)" ] }, { "cell_type": "markdown", "id": "08988978", "metadata": {}, "source": [ "We could use the standard approach based on the `dune.fem.operator.galerkin` construction\n", "but this requires adding the inverse mass operator later on. `dune-fem` provide an alternative approach\n", "where the operator automatically includes the inverse mass matrix:" ] }, { "cell_type": "code", "execution_count": 7, "id": "5c18dce7", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:35:01.269246Z", "iopub.status.busy": "2024-02-29T12:35:01.268850Z", "iopub.status.idle": "2024-02-29T12:35:45.469704Z", "shell.execute_reply": "2024-02-29T12:35:45.468293Z" } }, "outputs": [], "source": [ "from dune.fem.operator import molGalerkin as molOperator\n", "op = molOperator(form)" ] }, { "cell_type": "markdown", "id": "1473f80a", "metadata": {}, "source": [ "So `op(uh,res)` applies the given operator to `uh`,\n", "applies the inverse mass matrix and stores the result in `res`.\n", "Consequently in our case, `op(uh,res)` computes\n", "$M^{-1}\\vec{F}(\\vec{u}^n)$ as required.\n", "\n", "The remaining code is a rather standard time loop where we compute for a\n", "given old time step vector $\\vec{u}$ the new vector\n", "$$ \\frac{1}{2}\\Big( \\vec{u}_1^{n+1} + \\vec{u}_2^{n+1} \\Big) $$\n", "\n", "by applying the above operator twice.\n", "\n", "We show the transport of a a simple discontinuety:" ] }, { "cell_type": "code", "execution_count": 8, "id": "43386624", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:35:45.475646Z", "iopub.status.busy": "2024-02-29T12:35:45.475269Z", "iopub.status.idle": "2024-02-29T12:36:21.729561Z", "shell.execute_reply": "2024-02-29T12:36:21.728387Z" } }, "outputs": [ { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrAPADASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+imSmQQuYUR5QPlV2Kgn3IBx+RrBm8TtHpNhfJZBvtFi1/Ihmx5caqhYA7fmb5xgcZweRQBq39/9h+zgW01xJPL5SJDtznazZO5gMYU96h/tO7/6Aeof99wf/Haj1SeEavosBlQTNdOwjLDcQIJcnHXFTazrNnoOnNfXrMIgQoCLlmY9AP8A69NJydkKUlFXY3+07v8A6Aeof99wf/HaP7Tu/wDoB6h/33B/8drl/wDha+hf8+mo/wDftP8A4uj/AIWvoX/PpqP/AH7T/wCLrf6rW/lMPrdH+Y6j+07v/oB6h/33B/8AHaP7Tu/+gHqH/fcH/wAdrl/+Fr6F/wA+mo/9+0/+Lo/4WvoX/PpqP/ftP/i6Pqtb+UPrdH+Y6j+07v8A6Aeof99wf/HaP7Tu/wDoB6h/33B/8drl/wDha+hf8+mo/wDftP8A4uj/AIWvoX/PpqP/AH7T/wCLo+q1v5Q+t0f5jqP7Tu/+gHqH/fcH/wAdo/tO7/6Aeof99wf/AB2naNrNnr2nLfWTMYiSpDrhlYdQf/rVoVg04uzN4yUldGb/AGnd/wDQD1D/AL7g/wDjtH9p3f8A0A9Q/wC+4P8A47WlRSGZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1pUUAZv9p3f/QD1D/vuD/47UllqJu7qe2ksri1lhRJCJih3KxYDGxm/uHrV6s2D/kZb7/rzt//AEOagC/LFHPE8UsayRuCrI4yGB6gjvVSPRNKigjgj0yySKNxIka26hVYAAMBjg4AGfar1FAGbqf/AB/6N/1+N/6IlqxqGnWmq2T2l7Cs0D43ISR09xyKr6n/AMf+jf8AX43/AKIlrSpptO6E0mrM5j/hXvhb/oF/+TEv/wAVR/wr3wt/0C//ACYl/wDiq6eitPb1f5n95n9Xpfyr7kcx/wAK98Lf9Av/AMmJf/iqP+Fe+Fv+gX/5MS//ABVdPRR7er/M/vD6vS/lX3I5j/hXvhb/AKBf/kxL/wDFUH4e+Fsf8gv/AMmJf/iq6eg9KPb1f5n94fV6X8q+5GP4Ut4bbwppawxhA1rHIwHdmUEn8Sa2Ky/Df/IraR/15Q/+gCtSpqfGyqfwL0CiiioLCiiigAooooAKKKKACiiigAooooAKKKKACiiigArNg/5GW+/687f/ANDmrSrNg/5GW+/687f/ANDmoA0qKKKAM3U/+P8A0b/r8b/0RLWlWbqf/H/o3/X43/oiWtKgAooooAKKKKACg9KKD0oAy/Df/IraR/15Q/8AoArUrL8N/wDIraR/15Q/+gCtSrqfGyKfwL0CiiioLCiiigAooooAKKKKACiiigAooooAKKKKACiiigArNg/5GW+/687f/wBDmrSrNg/5GW+/687f/wBDmoAl1e8ksNKuLqIKXjXILjKrzgs3sOp6cCuYvfGNzBp1pNC9lLI088crAZVhG2FKguMblKsAC7YYbVfrXaUUAZOqTKNX0WHD7zdOwIjbbjyJf4sYB9s5rWrN1P8A4/8ARv8Ar8b/ANES1pUAFFFFABRRRQAUHpRQelAGX4b/AORW0j/ryh/9AFalZfhv/kVtI/68of8A0AVqVdT42RT+BegUUUVBYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFZsH/Iy33/AF52/wD6HNWlWbB/yMt9/wBedv8A+hzUAaVFFFAGbqf/AB/6N/1+N/6IlrSrN1P/AI/9G/6/G/8AREtaVABRRRQAUUUUAFB6UUHpQBl+G/8AkVtI/wCvKH/0AVqVl+G/+RW0j/ryh/8AQBWpV1PjZFP4F6BRRRUFhRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAVmwf8AIy33/Xnb/wDoc1aVZsH/ACMt9/152/8A6HNQBpUVS1eS6i0q4ezDeeF+Uou5gM8lV5yQMkDByR0Nc3calr/2G1NuLvcks4kZrI+ZNtceShG3CB0PL4AUjnByKAN/U/8Aj/0b/r8b/wBES1pVk6pIw1fRY/JcobpyZQV2g+RLx1zn8MVrUAFFFFABRRRQAUHpRQelAGX4b/5FbSP+vKH/ANAFalZfhv8A5FbSP+vKH/0AVqVdT42RT+BegUUUVBYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFZsH/Iy33/Xnb/8Aoc1aVZsH/Iy33/Xnb/8Aoc1AGlRRRQBm6n/x/wCjf9fjf+iJa0qzdT/4/wDRv+vxv/REtaVABRRRQAUUUUAFB6UUHpQBl+G/+RW0j/ryh/8AQBWpWX4b/wCRW0j/AK8of/QBWpV1PjZFP4F6BRRRUFhRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAVmwf8jLff9edv/6HNWlWbB/yMt9/152//oc1AGlRRRQBm6n/AMf+jf8AX43/AKIlrSrN1P8A4/8ARv8Ar8b/ANES1pUAFFFFABRRRQAUHpRQelAGX4b/AORW0j/ryh/9AFalZfhv/kVtI/68of8A0AVqVdT42RT+BegUUUVBYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFZsH/Iy33/AF52/wD6HNWlWbB/yMt9/wBedv8A+hzUAS6vHdS6VcR2ZYTlflCNtYjPIVuMEjIByME9RXPz2viSTS7OK2jkjaG+8yRZrzbI8ImVkXeu7cojJDAnJK/xA/N1tFAGTqhm/tfRQETyftT5cudwPkS4G3GMe+fwrWrN1P8A4/8ARv8Ar8b/ANES1pUAFFFFABRRRQAUHpRQelAGX4b/AORW0j/ryh/9AFalZfhv/kVtI/68of8A0AVqVdT42RT+BegUUUVBYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFZsH/Iy33/AF52/wD6HNWlWbB/yMt9/wBedv8A+hzUAaVFFFAGbqf/AB/6N/1+N/6IlrSrN1P/AI/9G/6/G/8AREtaVABRRRQAUUUUAFB6UUHpQBl+G/8AkVtI/wCvKH/0AVqVl+G/+RW0j/ryh/8AQBWpV1PjZFP4F6BRRRUFhRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAVmwf8AIy33/Xnb/wDoc1aVZsH/ACMt9/152/8A6HNQBpUVS1ezkv8ASri1iKh5FwA5wrc5Kt7HofY1gt4fv3sLKzay0wpDePMW81t0UZfzFWI+X8pBO3jHyqMYJ+UA2dT/AOP/AEb/AK/G/wDREtaVZOqRsdX0WTznCC6cGIBdpPkS89M5/HFa1ABRRRQAUUUUAFB6UUHpQBl+G/8AkVtI/wCvKH/0AVqVl+G/+RW0j/ryh/8AQBWpV1PjZFP4F6BRRRUFhRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAVmwf8AIy33/Xnb/wDoc1aVZsH/ACMt9/152/8A6HNQBpUUUUAZup/8f+jf9fjf+iJa0qzdT/4/9G/6/G/9ES1pUAFFFFABRRRQAUHpRQelAGX4b/5FbSP+vKH/ANAFalZfhv8A5FbSP+vKH/0AVqVdT42RT+BegUUUVBYUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFZsH/Iy33/Xnb/8Aoc1aVZsH/Iy33/Xnb/8Aoc1AGlRTJZY4InlmkWONFLO7nAUDqST0FQNqenrDDM19bCKdgsLmVdshPQKc8n6UAV9T/wCP/Rv+vxv/AERLWlWbqf8Ax/6N/wBfjf8AoiWtKgAooooAKKKKACg9KKbI6xxs7sFRQSzMcAAdzQBm+G/+RW0j/ryh/wDQBWpXNeH9f0aHw3pcUur2EciWkSsjXKAqQgyCM8GtL/hJNC/6DWnf+BSf41rOEuZ6GVOpHkWpp0Vmf8JJoX/Qa07/AMCk/wAaP+Ek0L/oNad/4FJ/jU8kuxXtIdzTorM/4STQv+g1p3/gUn+NH/CSaF/0GtO/8Ck/xo5Jdg9pDuadFZn/AAkmhf8AQa07/wACk/xo/wCEk0L/AKDWnf8AgUn+NHJLsHtIdzTorM/4STQv+g1p3/gUn+NH/CSaF/0GtO/8Ck/xo5Jdg9pDuadFZn/CSaF/0GtO/wDApP8AGj/hJNC/6DWnf+BSf40ckuwe0h3NOisz/hJNC/6DWnf+BSf40f8ACSaF/wBBrTv/AAKT/Gjkl2D2kO5p0Vmf8JJoX/Qa07/wKT/Gj/hJNC/6DWnf+BSf40ckuwe0h3NOisz/AISTQv8AoNad/wCBSf40f8JJoX/Qa07/AMCk/wAaOSXYPaQ7mnWbB/yMt9/152//AKHNVq0v7O/jaSzu4LlFOC0MgcA+mRVWD/kZb7/rzt//AEOapatuUmnqifUrL+0NPltfM8svgq+3IDAgjI7jIGR3rGk8LzzWphfUI8uLlJGW3IzHcMHkAG/htwO1ucA4IPU9JRSGZOqQqdX0WbL7xdOoAkbbjyJf4c4J98ZrWqjqNlPdtaSW1xHDLbTGUGSIyK2UdMYDL/fz17VXni19beVoL3TpJQhKI1m4DNjgE+bxzQBrUVm+Rrf/AEENP/8AAF//AI7R5Gt/9BDT/wDwBf8A+O0AaVFZvka3/wBBDT//AABf/wCO0eRrf/QQ0/8A8AX/APjtAGlWfrmlR65otzpskjRrOoG8DO0ggg478gcVDNFr6oDDe6c7b1BBs3HylhuP+t7DJ/CpPI1v/oIaf/4Av/8AHaabi7oUkpKzOE/4VD/1HP8AyU/+zo/4VD/1HP8AyU/+zru/I1v/AKCGn/8AgC//AMdo8jW/+ghp/wD4Av8A/Ha6frtf+b8Ecv1HD/y/i/8AM4T/AIVD/wBRz/yU/wDs6P8AhUP/AFHP/JT/AOzru/I1v/oIaf8A+AL/APx2o5ItfDxCO905lL4kJs3G1dp5H73nnA/Gj67X/m/BB9Rw/wDL+L/zOI/4VD/1HP8AyU/+zo/4VD/1HP8AyU/+zru/I1v/AKCGn/8AgC//AMdo8jW/+ghp/wD4Av8A/HaPrtf+b8EH1HD/AMv4v/M4T/hUP/Uc/wDJT/7Oj/hUP/Uc/wDJT/7Ou78jW/8AoIaf/wCAL/8Ax2jyNb/6CGn/APgC/wD8do+u1/5vwQfUcP8Ay/i/8zhP+FQ/9Rz/AMlP/s6P+FQ/9Rz/AMlP/s67cxa+LhFF7pxiKMWf7G+Q2RgY83uC35e9SeRrf/QQ0/8A8AX/APjtH12v/N+CD6jh/wCX8X/mcJ/wqH/qOf8Akp/9nR/wqH/qOf8Akp/9nXd+Rrf/AEENP/8AAF//AI7R5Gt/9BDT/wDwBf8A+O0fXa/834IPqOH/AJfxf+Zwn/Cof+o5/wCSn/2dH/Cof+o5/wCSn/2dd35Gt/8AQQ0//wAAX/8AjtRiLXzcOpvdOEQRSr/Y3yWycjHm9gF/P2o+u1/5vwQfUcP/AC/i/wDM4j/hUP8A1HP/ACU/+zo/4VD/ANRz/wAlP/s67vyNb/6CGn/+AL//AB2jyNb/AOghp/8A4Av/APHaPrtf+b8EH1HD/wAv4v8AzOE/4VD/ANRz/wAlP/s6P+FQ/wDUc/8AJT/7Ou78jW/+ghp//gC//wAdo8jW/wDoIaf/AOAL/wDx2j67X/m/BB9Rw/8AL+L/AMzH8JeCo/C01xP9ta6lmUJny9gUZz0yc1sQf8jLff8AXnb/APoc1Rxxa+XlEl7pyqHxGRZudy7Ryf3vHOR+FTWVldw39xd3d1DM8sUcQEMBjChC57u2c7/bpWE5yqS5pPU6KdONOPLFaGhRRRUFhWBc6Ndt4nh1OJbZ0WRWLu5WQJsZGjGFPy5bf15YYx3G/RQBheHdGuNI87zvIG6OOM+SSfNdd26Z8gfO+4Z6/dHJrdoooAKKKKAObTRL2DW7q9WCynjmSdW8yRlaYPtKhxsIwuzYOT8rE9sHQ8PafNpmkrb3EVvFMXZ3S2YmIFjk7RtXaOeBj6knJOpRQAUUUUAFcrZ+HdStlv4UmtbcXUaoZ4yzM7LIxLsuF+Z1cgkNkbRgnt1VFAGdoFhNpXh7T9PuGhaW2t0hYwqVT5QBwD24/wD1Vo0UUAFRXUTT2k0KSGJ5EZVkXqpIxkfSpaKAON/4RW9fSnsZLeyW3Nyk/wBkiuHSNwIthQkJlRuAkyAckke56uxhlt7C2hnkEs0cSo8gXaGYAAkDtk9qnooAKKKKAM3W9Pk1GyjiSOCYJMsjQXBxHKB/Cxwfr0PIFYS+GdSa80SaV7VjYQxRO/mEk7GOWGUzllxwGTnIYuOK6+igAooooAKy/Eeny6r4evrGCG3lmnhZIxcMVRWIwGyFY5HUcdRWpRQBz13o99farbX0sdojJ5XIlZmt9khZvLOwZ8xSFb7uAP4q6GiigD//2Q==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrANoDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+iiigCnf3/2H7OBbTXEk8vlIkO3OdrNk7mAxhT3qH+07v/oB6h/33B/8do1P/j/0b/r8b/0RLWlQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlRQBm/2nd/9APUP++4P/jtH9p3f/QD1D/vuD/47WlWLrHirSdBu7e1v52SSbkYQkKvPzH2yMVUYyk7RVyZSjFXk7Fj+07v/AKAeof8AfcH/AMdo/tO7/wCgHqH/AH3B/wDHa0UdZEV0YMrDIYHII9aWpKM3+07v/oB6h/33B/8AHaP7Tu/+gHqH/fcH/wAdrSooAzf7Tu/+gHqH/fcH/wAdo/tO7/6Aeof99wf/AB2tKigDNXVpPtdtBPpd5b/aHMaSSNEVDBGbna5PRT2rSrN1P/j/ANG/6/G/9ES1pUAQ3UnlWksgmig2oT5swyie7cjj8RXLN4vkk0LTL62lsXkntDNKASytMERvs6Ybh23nGcnjoa6+igDJ1SZRq+iw4febp2BEbbceRL/FjAPtnNVvFfiiPwvZQzvaSXLzOVRVbavHJy2Dj8qu6n/x/wCjf9fjf+iJauXVna30Pk3dtDcRZzslQOM/Q1UHFSTkroiak4tRdmecf8Le/wCoH/5N/wD2FH/C3v8AqB/+Tf8A9hXdf8I3oX/QF07/AMBU/wAKP+Eb0L/oC6d/4Cp/hXX7XDfyficnssV/z8/A4X/hb3/UD/8AJv8A+wo/4W9/1A//ACb/APsK7r/hG9C/6Aunf+Aqf4Uf8I3oX/QF07/wFT/Cj2uG/k/EPZYr/n5+Bwv/AAt7/qB/+Tf/ANhR/wALe/6gf/k3/wDYV3X/AAjehf8AQF07/wABU/wo/wCEb0L/AKAunf8AgKn+FHtcN/J+IeyxX/Pz8Dhf+Fvf9QP/AMm//sKP+Fvf9QP/AMm//sK7r/hG9C/6Aunf+Aqf4Uf8I3oX/QF07/wFT/Cj2uG/k/EPZYr/AJ+fgcL/AMLe/wCoH/5N/wD2FH/C3v8AqB/+Tf8A9hXdf8I3oX/QF07/AMBU/wAKP+Eb0L/oC6d/4Cp/hR7XDfyfiHssV/z8/A4X/hb3/UD/APJv/wCwo/4W9/1A/wDyb/8AsK7r/hG9C/6Aunf+Aqf4Uf8ACN6F/wBAXTv/AAFT/Cj2uG/k/EPZYr/n5+Bwv/C3v+oH/wCTf/2FH/C3v+oH/wCTf/2Fd1/wjehf9AXTv/AVP8KP+Eb0L/oC6d/4Cp/hR7XDfyfiHssV/wA/PwOF/wCFvf8AUD/8m/8A7CuP1rVJ/F/iWOWOFonuDHBFC8u8IeAADgYBJz9Sa9q/4RvQv+gLp3/gKn+FYniuPSvDegveW2mW0UpnhC+RCqFisivgkDp8laUa9GM/3cNXpuZVqFaUP3s9FrsdPZwtbWNvAxUtHGqEqMDIGOB2qeqWk6lFrGlW+oQK6xzLuCv1HOP6Vdrz5JptM9KLTSa2CiiikMKKKKAM3U/+P/Rv+vxv/REtaVZup/8AH/o3/X43/oiWtKgAooooAzdT/wCP/Rv+vxv/AERLWlWbqf8Ax/6N/wBfjf8AoiWtKgAooooAKKKKACiiigAooooAKKKKACiiigAooooAKwvGNjBf+F7yO4XKxqJFI6gg9f5/nW7WX4l/5FjVD6Wsh/JSaum7TTXciqk4NPsW9PsYNMsILK1XbDCu1Qev4+9WaKKltt3ZSSSsgooopDCiiigDN1P/AI/9G/6/G/8AREtaVZup/wDH/o3/AF+N/wCiJa0qAIrm5hs7aS4nfZFGMscE/kByT7DrVGXxBpkNtDcPO4jlLhT5LkrsOHLDblQp4JbGO9W76zjv7OS2lZlVwPmQ/MpByCM9wQD+FY1x4Qs7u3ihnu7uQRTSTKzCMkM53MR8nyndkhhhl3HBA4oAv6n/AMf+jf8AX43/AKIlrSrJ1SFTq+izZfeLp1AEjbceRL/DnBPvjNa1ABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAVm+IV3eGtVX1s5h/44a0qx/Fcs8HhXU3tohJJ5DDaf7p4Y/gCT+FXT1miKmkH6GtE26JG9VBp1c94K1HUNU8Nw3GoxBJQxRGC7d6jADY+ufyroaU4uMnF9AhJTipLqFFFFSWFFFFAGbqf/H/o3/X43/oiWtKs3U/+P/Rv+vxv/REtaVABRRRQBm6n/wAf+jf9fjf+iJa0qzdT/wCP/Rv+vxv/AERLWlQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFVdSG7S7sesLj/x01aqC8G6xuB6xsP0prcUtmR6Wd2k2R9YEP8A46Kt1R0Y7tC08+ttGf8Ax0Vepy+Jij8KCiiipKCiiigDN1P/AI/9G/6/G/8AREtaVZup/wDH/o3/AF+N/wCiJa0qACiqWrx3UulXEdmWE5X5QjbWIzyFbjBIyAcjBPUVzdxZ+IZbG1SEX0ZhlnIQ3K7my4aDzG35aMKSrDJYkcZ4NAG/qf8Ax/6N/wBfjf8AoiWtKsnVDN/a+igInk/any5c7gfIlwNuMY98/hWtQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFMmG6CQeqkfpT6QjII9aAM/QDu8OaWfW0iP8A44K0ay/DZz4X0n/rzh/9AFalVP4mRT+BBRRRUlhRRRQBm6n/AMf+jf8AX43/AKIlrSrN1P8A4/8ARv8Ar8b/ANES1pUAFFFFAGbqf/H/AKN/1+N/6IlrSrN1P/j/ANG/6/G/9ES1pUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABTZJEijaSR1REBZmY4AA6kmnVHPBHc28kEyh4pUKOp7gjBFCBmJ4P1Kzv/DdnHa3CSvbwpFKo6owGOR+H41v1zPgjRLPSdDWW2DmS6w8rOcnIyAB7Dn8zXTVpW5faPl2MqPN7Nc24UUUVmahRRRQBm6n/wAf+jf9fjf+iJa0qzdT/wCP/Rv+vxv/AERLWlQAUVU1K9/s/T5boR+YyYCpu2gsSAMnsMkZPasa88Xw2FjZy3EUST3F4bUxtOFVds3lO4YgZAJBAxk5HTkgA0tT/wCP/Rv+vxv/AERLWlWTqk8I1fRYDKgma6dhGWG4gQS5OOuK1qACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAy/Dn/IBth6bh+TGtSsvw7/AMgSIekko/KRq1Kup8bIp/AvQKKKKgsKKKKAM3U/+P8A0b/r8b/0RLWlWbqf/H/o3/X43/oiWtKgBksUc8TxSxrJG6lXRxkMD1BHcVDFp1lBAYYbO3jiLq5RIlC7hjBwB1G1cH2HpVmigDN1P/j/ANG/6/G/9ES1pVm6n/x/6N/1+N/6IlrSoAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigDL8Pf8gjHpcXA/KZ61Ky/D//ACDZR6Xl0P8AyPJWpV1PjZFP4F6BRRRUFhRRRQBm6n/x/wCjf9fjf+iJa0qzdT/4/wDRv+vxv/REtaVABRTJZY4InlmkWONFLO7nAUDqST0FV21XTljgka/tQlwdsDGZcSnOMKc88kdKAINT/wCP/Rv+vxv/AERLWlWbqf8Ax/6N/wBfjf8AoiWtKgAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKAOQ8G+JbbVLrUNOihkR4riadWbGGRpSfwPzdK6+uZ8K6HY6fc6neW8ZWaW6liPPCqHOAPTt+Qrpq1r8vO+XYxoc/s1z7hRRRWRsFFFFAGbqf/H/o3/X43/oiWtKs3U/+P/Rv+vxv/REtaVAFTUrL+0NPltfM8svgq+3IDAgjI7jIGR3rKOg35tIrU6lbmAXb3M6fZG/e7n8zZ/rOBvLH3G0HODu6CigDJ1SNjq+iyec4QXTgxALtJ8iXnpnP44rWrN1P/j/0b/r8b/0RLWlQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFAGXon3dQHpfS/zB/rWpWXov+s1Qel8/wD6Cp/rWpVz+Iin8IUUUVBYUUUUAZup/wDH/o3/AF+N/wCiJa0qzdT/AOP/AEb/AK/G/wDREtaVABRRRQBm6n/x/wCjf9fjf+iJa0qzdT/4/wDRv+vxv/REtaVABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAZej/8fOrj0vj/AOi4z/WtSvPfAOq3994i1uK6uHkj3mQq0YX5shc+xwoGPb2r0Ktq8HCdn5GNCanDmXmFFFFYmwUUUUAZup/8f+jf9fjf+iJa0qzdT/4/9G/6/G/9ES1pUAUtXkuotKuHsw3nhflKLuYDPJVeckDJAwckdDXOzajrZ06Iw/bvNU3IVjZ4aVww8hXBT5UZDlmAXBGMqciuvooAydUkYavosfkuUN05MoK7QfIl465z+GK1qzdT/wCP/Rv+vxv/AERLWlQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFAHmHjeK/D39ppluSt1eIZhAp3v+4BC4HUcOT7iub0vx94h0uKKNblbi3jBVY7hNw/76GG4+tbOuTXeufElbLSriW3kikCli23DIMMw/DP1r1qvSlVjSpxjOKd0eXGlKrUlKEmrOx49/wtfXf+fTTv8Av2//AMXR/wALX13/AJ9NO/79v/8AF17DRWX1ij/z7/H/AIBt9Wrf8/fw/wCCePf8LX13/n007/v2/wD8XWv4Z+Imr6xr1vYT2Fq8cxwTArKyD+9ySMCvSqKmVek00qf4lRw9ZNN1PwM3U/8Aj/0b/r8b/wBES1pVm6n/AMf+jf8AX43/AKIlrSrkOwKKKKAM3U/+P/Rv+vxv/REtaVUdRsp7trSS2uI4ZbaYygyRGRWyjpjAZf7+evaq80WvqgMN7pztvUEGzcfKWG4/63sMn8KANais3yNb/wCghp//AIAv/wDHaPI1v/oIaf8A+AL/APx2gDSorN8jW/8AoIaf/wCAL/8Ax2jyNb/6CGn/APgC/wD8doA0qKyYYtfZCZr3TkbewAFm5+UMdp/1vcYP41J5Gt/9BDT/APwBf/47QBpUVm+Rrf8A0ENP/wDAF/8A47R5Gt/9BDT/APwBf/47QBpUVkmLXxcIovdOMRRiz/Y3yGyMDHm9wW/L3qTyNb/6CGn/APgC/wD8doA0qKzfI1v/AKCGn/8AgC//AMdo8jW/+ghp/wD4Av8A/HaANKismeLX1t5WgvdOklCEojWbgM2OAT5vHNSeRrf/AEENP/8AAF//AI7QBpUVm+Rrf/QQ0/8A8AX/APjtHka3/wBBDT//AABf/wCO0AY0/gSzm8XR68ty8e2QTNAF4Zx3znjnBIx6+tdXWb5Gt/8AQQ0//wAAX/8AjtRwxa+yEzXunI29gALNz8oY7T/re4wfxq51JTtzPYiFOML8q3Nais3yNb/6CGn/APgC/wD8do8jW/8AoIaf/wCAL/8Ax2oLNKis3yNb/wCghp//AIAv/wDHajMWvi4RRe6cYijFn+xvkNkYGPN7gt+XvQBJqf8Ax/6N/wBfjf8AoiWtKsv7BqM17ZzXd9avHbSmXZFashYlGTqZDj7+enatSgAooooAK5eXw9fnU7+aJ7eIXUNxEblXYS/vAuwkbedhXaPm6MTx0PUUUAZeh6dJpttOjxW9uss3mJb2xJjhG1VwpwvUqW6Dlj9a1KKKACsLxNo1xrNtHHB5BKrIuJiQEZlwsq4B+ZTyOnU8it2igDDh027fxMupz2dhEi2+zfBKTKzkDduOwblGAF59yMkBdyiigAqK6iae0mhSQxPIjKsi9VJGMj6VLRQBxk3hXU5/DOpaRG1hZpdfMsce6SMjyghQjCkAsquSM53MCDznsIg6xIsjKzhQGKjAJ74HOKfRQAUUUUAcvJ4f1A6nfzwyW8QuoriL7UrsJv3gXYxG3nyymAN3RieOh0fDulyaRp72zpDEhlLxwxOXWMEDI3EAtlgzcjjdjoK16KACiiigDC8TaNcazbRxweQSqyLiYkBGZcLKuAfmU8jp1PIqvFoV8fFyazKLZUKAOqPuYfu9pUHYCfm5zuC4/g3fNXS0UAFFFFABXMHRL8x6xbi2sIra8u1uUVJmw4URAo6+WAA+x9x5+90bmunooAo6PZSafpkdtLsDKzsEjOUjDOWCLwOFBCjgcDoOlXqKKAP/2Q==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "uh = space.interpolate(conditional(x[0]<0.25,1,0), name=\"solution\")\n", "uh.plot()\n", "w = uh.copy()\n", "un = uh.copy()\n", "cfl = 0.1\n", "tau = cfl/(speed*N)\n", "\n", "t = 0\n", "while t<1:\n", " un.assign(uh)\n", "\n", " op(uh, w)\n", " uh.axpy(-tau,w) # with numpy backend equivalent to 'uh.as_numpy[:] -= tau*w.as_numpy[:]'\n", " op(uh, w)\n", " uh.axpy(-tau,w) # with numpy backend equivalent to 'uh.as_numpy[:] -= tau*w.as_numpy[:]'\n", "\n", " # with numpy backend the following is equivalent to\n", " # 'uh.as_numpy[:] = 0.5*(uh.as_numpy[:] + un.as_numpy[:])'\n", " uh *= 0.5\n", " uh.axpy(0.5,un)\n", " t += tau\n", "uh.plot()" ] }, { "cell_type": "markdown", "id": "88ced265", "metadata": {}, "source": [ "\n", ".. todo:: add some remarks on spectral dg\n", "\n", "The `dune.femdg` module provides further methods for solving time\n", "dependent advection-diffusion problem. The module includes\n", "implicit-explicit Runge-Kutta methods and stabilization approaches for\n", "non-linear advection dominated models. This module is described in\n", "more details in [other places in this tutorial](femdg.rst)." ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }