diff --git a/pyro/multigrid/axisymmetric_MG.py b/pyro/multigrid/axisymmetric_MG.py new file mode 100644 index 000000000..c17523e5b --- /dev/null +++ b/pyro/multigrid/axisymmetric_MG.py @@ -0,0 +1,135 @@ +r""" +A multigrid solver for axisymmetric coordinates (r-z) to solve +an elliptic equation of the form: + +L \phi = f + +where L is the Laplacian operator in cylindrical coords. +""" + +import matplotlib.pyplot as plt +import numpy as np +import pyro.mesh.array_indexer as ai +from pyro.multigrid import MG + +np.set_printoptions(precision=3, linewidth=128) + + +class AxisymmetricMG2d(MG.CellCenterMG2d): + + def __init__(self, nx, ny, xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0, + xl_BC_type="reflect", xr_BC_type="dirichlet", + yl_BC_type="dirichlet", yr_BC_type="dirichlet", + xl_BC=None, xr_BC=None, + yl_BC=None, yr_BC=None, + nsmooth=10, nsmooth_bottom=50, + verbose=0, + true_function=None, vis=0, vis_title=""): + + # initialize the MG object with the auxillary "coeffs" field + MG.CellCenterMG2d.__init__(self, nx, ny, ng=1, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xl_BC_type=xl_BC_type, xr_BC_type=xr_BC_type, + yl_BC_type=yl_BC_type, yr_BC_type=yr_BC_type, + xl_BC=xl_BC, xr_BC=xr_BC, + yl_BC=yl_BC, yr_BC=yr_BC, + nsmooth=nsmooth, nsmooth_bottom=nsmooth_bottom, + verbose=verbose, + true_function=true_function, vis=vis, + vis_title=vis_title) + + def smooth(self, level, nsmooth): + """ + Use red-black Gauss-Seidel iterations + """ + + v = self.grids[level].get_var("v") + f = self.grids[level].get_var("f") + + self.grids[level].fill_BC("v") + + myg = self.grids[level].grid + + dr = myg.dx + dz = myg.dy + + rc = ai.ArrayIndexer(myg.x2d, grid=myg) + rl = rc - 0.5 * dr + rr = rc + 0.5 * dr + + denom = 2 * (1.0 / dr**2 + 1.0 / dz**2) + + # do red-black G-S + for _ in range(nsmooth): + + # do the red black updating in four decoupled groups + # + # + # | | | + # --+-------+-------+-- + # | | | + # | 4 | 3 | + # | | | + # --+-------+-------+-- + # | | | + # jlo | 1 | 2 | + # | | | + # --+-------+-------+-- + # | ilo | | + # + # groups 1 and 3 are done together, then we need to + # fill ghost cells, and then groups 2 and 4 + + for n, (ix, iy) in enumerate([(0, 0), (1, 1), (1, 0), (0, 1)]): + + v.ip_jp(ix, iy, s=2)[:, :] = (-f.ip_jp(ix, iy, s=2) + + 1.0 / (rc.ip(ix, s=2) * dr**2) * (rr.ip(ix, s=2) * v.ip_jp(1+ix, iy, s=2) + + rl.ip(ix, s=2) * v.ip_jp(-1+ix, iy, s=2)) + + 1.0 / dz**2 * (v.ip_jp(ix, 1+iy, s=2) + v.ip_jp(ix, -1+iy, s=2))) / denom + + if n in (1, 3): + self.grids[level].fill_BC("v") + + if self.vis == 1: + plt.clf() + + plt.subplot(221) + self._draw_solution() + + plt.subplot(222) + self._draw_V() + + plt.subplot(223) + self._draw_main_solution() + + plt.subplot(224) + self._draw_main_error() + + plt.suptitle(self.vis_title, fontsize=18) + + plt.draw() + plt.savefig("mg_%4.4d.png" % (self.frame)) + self.frame += 1 + + def _compute_residual(self, level): + """ compute the residual and store it in the r variable""" + + v = self.grids[level].get_var("v") + f = self.grids[level].get_var("f") + r = self.grids[level].get_var("r") + + myg = self.grids[level].grid + + dr = myg.dx + dz = myg.dy + + rc = ai.ArrayIndexer(myg.x2d, grid=myg) + rl = rc - 0.5 * dr + rr = rc + 0.5 * dr + + # compute the residual + # r = f - L_eta phi + L_phi = (1.0 / (rc.v() * dr**2) * (rr.v() * v.ip(1) - 2.0 * rc.v() * v.v() + rl.v() * v.ip(-1)) + + 1.0 / (dz**2) * (v.jp(1) - 2.0 * v.v() + v.jp(-1))) + + r.v()[:, :] = f.v() - L_phi diff --git a/pyro/multigrid/multigrid-axisymmetric.ipynb b/pyro/multigrid/multigrid-axisymmetric.ipynb new file mode 100644 index 000000000..f2b43e8ad --- /dev/null +++ b/pyro/multigrid/multigrid-axisymmetric.ipynb @@ -0,0 +1,880 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d276064d-6f7d-4a68-9c46-c0a41eef0646", + "metadata": {}, + "source": [ + "# Axisymmetric Multigrid Solver\n", + "\n", + "In axisymmetric coordinates ($r$-$z$), the Poisson equation appears as:\n", + "\n", + "$$\\frac{1}{r} \\frac{\\partial}{\\partial r} \\left ( r \\frac{\\partial \\phi}{\\partial r} \\right ) + \\frac{\\partial^2\\phi}{\\partial z^2} = f$$\n", + "\n", + "We will solve this with inhomogeneous Dirichlet boundary conditions on $[0, 1]^2$. With\n", + "\n", + "$$f = 4z^2 + 2 r^2$$\n", + "\n", + "and\n", + "\n", + "\\begin{align}\n", + "\\phi(r=0, z) &= 0 \\\\\n", + "\\phi(r=1, z) &= z^2 \\\\\n", + "\\phi(r, z=0) &= 0 \\\\\n", + "\\phi(r, z=1) &= r^2\n", + "\\end{align}\n", + "\n", + "the solution is:\n", + "\n", + "$$\\phi(r, z) = r^2 z^2$$" + ] + }, + { + "cell_type": "markdown", + "id": "9d8bd6c3-31e7-4eec-900f-0b9d09262316", + "metadata": {}, + "source": [ + "## Setting up the solver\n", + "\n", + "We'll take the convention that `x` $\\rightarrow$ `r`, and `y` $\\rightarrow$ `z` when accessing the grid data." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "2bd555f1-d90c-4ac5-862b-4699d7e2eee9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def f(r, z):\n", + " return 4 * z**2 + 2 * r**2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "af11c9e2-9c70-4075-b04a-c942e738bbf8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def rr_func(z):\n", + " return z**2\n", + "\n", + "def zr_func(r):\n", + " return r**2" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b796efb7-c279-449f-9862-0748d6940bb1", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def true(r, z):\n", + " return r**2 * z**2" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "8a663591-cc50-4042-9612-3514bbf472e2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import pyro.multigrid.axisymmetric_MG as MG" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "943a0e92-e969-4888-b5bf-827eb0f2251a", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "cc data: nx = 2, ny = 2, ng = 1\n", + " nvars = 3\n", + " variables:\n", + " v: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " f: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " r: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + "\n", + "cc data: nx = 4, ny = 4, ng = 1\n", + " nvars = 3\n", + " variables:\n", + " v: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " f: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " r: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + "\n", + "cc data: nx = 8, ny = 8, ng = 1\n", + " nvars = 3\n", + " variables:\n", + " v: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " f: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " r: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + "\n", + "cc data: nx = 16, ny = 16, ng = 1\n", + " nvars = 3\n", + " variables:\n", + " v: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " f: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " r: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + "\n", + "cc data: nx = 32, ny = 32, ng = 1\n", + " nvars = 3\n", + " variables:\n", + " v: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " f: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + " r: min: 0.0000000000 max: 0.0000000000\n", + " BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet \n", + "\n" + ] + } + ], + "source": [ + "N = 32\n", + "\n", + "mg = MG.AxisymmetricMG2d(N, N,\n", + " xl_BC_type=\"dirichlet\", xr_BC_type=\"dirichlet\",\n", + " yl_BC_type=\"dirichlet\", yr_BC_type=\"dirichlet\",\n", + " xr_BC=rr_func, yr_BC=zr_func, true_function=true,\n", + " verbose=1)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "77e19758-be49-4071-a8f8-cbdb2620e01e", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Source norm = 2.4028431443110723\n" + ] + } + ], + "source": [ + "mg.init_zeros()\n", + "mg.init_RHS(f(mg.x2d, mg.y2d))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "d6325d72-7261-4cba-a078-b9085c1f2ff8", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "source norm = 2.4028431443110723\n", + "<<< beginning V-cycle (cycle 1) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 2.4028431443110723\n", + " after G-S, residual L2: 12.949907537944354\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 8.851078613176675\n", + " after G-S, residual L2: 2.510200336301418\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 1.709803417983295\n", + " after G-S, residual L2: 0.3530965909891356\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 0.23964419976745627\n", + " after G-S, residual L2: 0.0015435410120936467\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 0.0018088549717517367\n", + " after G-S, residual L2: 1.0026872855338396e-06\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 0.7258102075427889\n", + " after G-S, residual L2: 0.007609680921336424\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 3.004878183732335\n", + " after G-S, residual L2: 0.03277677801146209\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 12.050024833019929\n", + " after G-S, residual L2: 0.19063830318826297\n", + "\n", + "cycle 1: relative err = 1.0000000000002403, residual err = 0.07933863832918714\n", + "\n", + "<<< beginning V-cycle (cycle 2) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 0.19063830318826297\n", + " after G-S, residual L2: 0.0926805547652613\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 0.06425201981038277\n", + " after G-S, residual L2: 0.02407070918527181\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 0.016455182610518017\n", + " after G-S, residual L2: 0.0043049260130849365\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 0.0029271166829855887\n", + " after G-S, residual L2: 1.8327710204196338e-05\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 2.1477835153368836e-05\n", + " after G-S, residual L2: 1.1907237712900042e-08\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 0.00897020761244355\n", + " after G-S, residual L2: 9.653883983912586e-05\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 0.049300464710413945\n", + " after G-S, residual L2: 0.00044614027356186585\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 0.16999414440442662\n", + " after G-S, residual L2: 0.0019202131720583244\n", + "\n", + "cycle 2: relative err = 10.978062371317007, residual err = 0.0007991421231987557\n", + "\n", + "<<< beginning V-cycle (cycle 3) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 0.0019202131720583244\n", + " after G-S, residual L2: 0.0010219200619397246\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 0.0007140734781043038\n", + " after G-S, residual L2: 0.0003093534580562075\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 0.00021337284756819334\n", + " after G-S, residual L2: 5.7468382932972244e-05\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 3.9129258918991144e-05\n", + " after G-S, residual L2: 2.3425278619161356e-07\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 2.745133430982456e-07\n", + " after G-S, residual L2: 1.522144018885369e-10\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 0.00012089821538130831\n", + " after G-S, residual L2: 1.3383258289960475e-06\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 0.001037508115565216\n", + " after G-S, residual L2: 1.0607738314113014e-05\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 0.0052861829869032435\n", + " after G-S, residual L2: 4.9888437474242e-05\n", + "\n", + "cycle 3: relative err = 2.09458302545174, residual err = 2.0762253080213313e-05\n", + "\n", + "<<< beginning V-cycle (cycle 4) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 4.9888437474242e-05\n", + " after G-S, residual L2: 2.538456722323062e-05\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 1.7843180248288583e-05\n", + " after G-S, residual L2: 7.693029414602172e-06\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 5.35608506733803e-06\n", + " after G-S, residual L2: 1.1551104371197595e-06\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 7.87278670221068e-07\n", + " after G-S, residual L2: 4.4297763451188705e-09\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 5.191043898074973e-09\n", + " after G-S, residual L2: 2.8790041126396894e-12\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 2.43932192249724e-06\n", + " after G-S, residual L2: 2.7890564989383105e-08\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 2.9949913505872747e-05\n", + " after G-S, residual L2: 3.3530860413328087e-07\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 0.00018118330298256175\n", + " after G-S, residual L2: 1.7285820461903024e-06\n", + "\n", + "cycle 4: relative err = 0.24059945376945252, residual err = 7.193902982318517e-07\n", + "\n", + "<<< beginning V-cycle (cycle 5) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 1.7285820461903024e-06\n", + " after G-S, residual L2: 8.644757393660188e-07\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 6.081052158846044e-07\n", + " after G-S, residual L2: 2.5731805841061746e-07\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 1.7954267840324867e-07\n", + " after G-S, residual L2: 3.454727504621579e-08\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 2.35532834401123e-08\n", + " after G-S, residual L2: 1.2802396474748287e-10\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 1.5002367951056565e-10\n", + " after G-S, residual L2: 8.321546399995043e-14\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 7.29271394181447e-08\n", + " after G-S, residual L2: 8.487866427430119e-10\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 1.0173339643986725e-06\n", + " after G-S, residual L2: 1.1609265675171233e-08\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 6.390999799107682e-06\n", + " after G-S, residual L2: 6.14053659923441e-08\n", + "\n", + "cycle 5: relative err = 0.008331101466533196, residual err = 2.5555295250015103e-08\n", + "\n", + "<<< beginning V-cycle (cycle 6) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 6.14053659923441e-08\n", + " after G-S, residual L2: 3.0550123706116485e-08\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 2.149229667092813e-08\n", + " after G-S, residual L2: 9.144022302194667e-09\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 6.38497722264134e-09\n", + " after G-S, residual L2: 1.2303742476680692e-09\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 8.389582140305002e-10\n", + " after G-S, residual L2: 4.52744417393934e-12\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 5.305420254183442e-12\n", + " after G-S, residual L2: 2.9429312028616896e-15\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 2.598009323501818e-09\n", + " after G-S, residual L2: 3.041147207346802e-11\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 3.679146069734179e-08\n", + " after G-S, residual L2: 4.1823570088659867e-10\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 2.303549682963893e-07\n", + " after G-S, residual L2: 2.2113032601875124e-09\n", + "\n", + "cycle 6: relative err = 0.0003640724193854146, residual err = 9.202861474428548e-10\n", + "\n", + "<<< beginning V-cycle (cycle 7) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 2.2113032601875124e-09\n", + " after G-S, residual L2: 1.0985112981199566e-09\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 7.728639267937154e-10\n", + " after G-S, residual L2: 3.3290872814456365e-10\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 2.3257305304952835e-10\n", + " after G-S, residual L2: 4.6050544673307166e-11\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 3.140392361594173e-11\n", + " after G-S, residual L2: 1.6964280734776735e-13\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 1.9879316605375316e-13\n", + " after G-S, residual L2: 1.1027163641068267e-16\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 9.729587284099105e-11\n", + " after G-S, residual L2: 1.1405976655701914e-12\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 1.3626014840494129e-09\n", + " after G-S, residual L2: 1.538230110731059e-11\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 8.437345807443207e-09\n", + " after G-S, residual L2: 8.074879380280347e-11\n", + "\n", + "cycle 7: relative err = 1.4977300083815499e-05, residual err = 3.36055201913545e-11\n", + "\n", + "<<< beginning V-cycle (cycle 8) >>>\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 8.074879380280347e-11\n", + " after G-S, residual L2: 4.0107708458269646e-11\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 2.8219552226511976e-11\n", + " after G-S, residual L2: 1.2298073043642456e-11\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 8.594525845722229e-12\n", + " after G-S, residual L2: 1.7431864268113942e-12\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 1.1888381167532685e-12\n", + " after G-S, residual L2: 6.43470558359894e-15\n", + "\n", + " bottom solve:\n", + " level: 0, grid: 2 x 2\n", + "\n", + " level: 1, grid: 4 x 4\n", + " before G-S, residual L2: 7.54040012905408e-15\n", + " after G-S, residual L2: 4.182692295487174e-18\n", + "\n", + " level: 2, grid: 8 x 8\n", + " before G-S, residual L2: 3.684697011767843e-12\n", + " after G-S, residual L2: 4.321110684251904e-14\n", + "\n", + " level: 3, grid: 16 x 16\n", + " before G-S, residual L2: 5.09793376542744e-11\n", + " after G-S, residual L2: 5.723202162194525e-13\n", + "\n", + " level: 4, grid: 32 x 32\n", + " before G-S, residual L2: 3.125457048392291e-10\n", + " after G-S, residual L2: 2.9832702148195685e-12\n", + "\n", + "cycle 8: relative err = 5.979790854009595e-07, residual err = 1.2415584520707083e-12\n", + "\n" + ] + } + ], + "source": [ + "mg.solve()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "8e84724b-2c51-4d88-b243-e892547f445a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "v = mg.get_solution()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "f6ea0519-6c94-4074-95b8-09d8fd4fe42d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "b246b1be-a7be-4070-a56e-fa2d1242350c", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAD6CAYAAAAm/xuVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAaoElEQVR4nO3db4wd13nf8e/vLpekSElxDNmJLRGVDDhGCCO2DEJyqiKI/4Zmg/hNXkhG/CJIwDdxKgcOCrsF6sZAkQQojPiFUJSI1LSIbMO1RUAIWP15YUMwECgiFdoWTTmgFSWiqIZmXNW0kVC7d5++uHfV5e6c587OnXvvzO7vAwy4vHNn5uzs3bNnznnOcxQRmJn10WDRBTAza8oVmJn1liswM+stV2Bm1luuwMyst/bM+4J7tS/2c7C9E6rZTik5MDtn6bj0fMnfiey4QXlfpMcVjknOl33P2XGRfGulfU2OAdI/t/lx1SPtWiqPwA8Ga8V9y8m+vUurxX37VN63f7BS/Xp2TPK5OvPta1ci4k3FN0zwK+87GP/4w2Gt95759rXHI+Jo02tNY+4V2H4Ocrc+sP0DC7+wWloqH5Psy47TnuS27F2uPma5+nUA9u0t7orC+QBif7JvX3nfcH91+Yc3lL/n4b7yL8PKgfK+1RvKldvqDYXXD2z/GIDVg+UKZ/VAeV8crK4Elg9WVxoAN934T8V9P3Pjj4v7Dh38P8V9b7vhSnHfO/a/Uvn64b3/UDzm55bLDYGlt1z4u+LOGq78cMjTj99W673Lb/n+LdNcaxpzr8DMrA+CYZRbml3hCszMtghgje4HubsCM7NKa7gFZmY9FAQrfoQ0sz4KYOhHyNmLtfJNlpIfQDJ4mV+wcM5sUvww+Uu2Vt6n1fK+WEqOK1xvsFIe/VvbU963tFL+3taST9Bgqfqc2cBxpPuahXMMCxdcScIofrK0r7jvSiEsA2CQfOaWkn0DFX5m6WPcD5J90+tDH9jEQFZJ+yX9laRvSTon6Q827Dsq6XuSLkj69GyLambzEsAwota2SHVaYNeA90fEjyUtA9+U9L+AZ4AHgA8BF4FnJD0aEd+dXXHNbF663wNWowKLUcKw9ci95fEWwF3AhYh4AUDSl4GPAq7AzHouiF70gdWaCylpSdJZ4DLwZEQ8DdwKvLThbRfHr1Udf1zSaUmnV7g2ZZHNbNYiYKXmtki1OvEjYgi8W9IbgJOS3kn17LnKbyciTgAnAG7WG7tfrZvtemKYTzTuhG1lo4iIV4FvAEcZtbgObdh9G3CprYKZ2eIEsBb1tkWa2AKT9CZgJSJelXQD8EHgjxl14r9d0h3Ay8C9wMdmWdhWJeEL2ToBKvzE0mOyfUmIRST71CD8Qkvlv1eD1WQydPKcsJSEX5RCImJPEk7wWvl8g2bz9olB9fc9TE742qA8Wf5q0jBZSjJVZCEWTayl7Y/qyeHb0YcWWJ1HyLcA/13SEqMW21ci4i8AJH0CeJxRVNVDEXFuZiU1s7kZBbLugAosIr4N3FnYdwo41XahzGyxAlhJk6x1Q+8j8c2sfYEY9iBhsyswM6u0FjvgEdLMdp8d0we2KyUjfJRSjGTjydl8sWxfw4ne2lMYhUxGNdOJ3skk6kEyQlka5CtN8gZYSiZK5zn9s4nehXMWRicBhoPyr8a1pBj/V0lO7JatpX1UZ6c8uxi6D8zM+miUkdUVmJn1UIR4Lctt1BGuwMys0pr7wMysj0ad+H6ENLNecie+mfWUO/G7IFlVJZIgvUZP/tkKLsNkifZkKJ9k0nOTEAtlefSzid4rWW7+JCSiNJk7+ZazUImkiOk5lwrnTMuh8s615Nfmn8unbCQLJh3OONB01udvw86uwMyskUCsRPerh+6X0Mzmri+d+N0voZnNXSCGUW+ro84KZpJ+b7zy2XOSviRp/6TzugIzs0prDGptk4xzCT4AfAQ4DNwn6fCm99wK/BvgSES8k1GOwXsnnduPkGa2RQTbCaO4RdLpDf8/MV4HY13dFcz2ADdIWgEOUCNFvSswM9ti1IlfeyrRlYg4kuyvWsHs7uuuF/GypP8M/D3wT8ATEfHEpAu7AquSZogo7EuOmUm+/CT3OsPqc2bny7JbZNkjIg2xqP4LXgprgDy0gcaZKkqvNwwTmEGIRSmsZ7iWXGvWYRTt9TBNXMFM0k8zapXdAbwK/E9JvxERf56d2H1gZrZFINai3lZDnRXMPgj8bUT8ICJWgEeAfznpxK7AzKzSkEGtrYbXVzCTtJdR5/yjm97z98B7JR2QJOADwPlJJ3YFZmZbjNaFHNTaJp4rYhVYX8HsPKOVzc4BSDol6a0R8TTwVeBZ4DuM6qYThVO+zn1gZlah3ZW5SyuYRcSxDV9/Fvjsds7rCszMthgtq+aEhmbWQxGq9Xi4aK7AtqkUEqFsUY9sX5apIku/kGWjKJ1zNVmcIwnLyDJODNKQiOrvu7jIBpMyTiSZKpTc41K4RPqE1PDxKfmlz0Isrm0/Oofh2qyzUeyACkzSIeB/AD8LrDGKsv3CeN+LwFVgCKxOCGYzs54Y5QPbGel0VoFPRcSzkm4Czkh6MiLWpwG8LyKuzK6IZjZ/OyQja0S8Arwy/vqqpPOMpgZsnsdkZjvEKIxiZ7TAXifpduBO4OnxSwE8ISmA/7ppAufG444DxwH2c6BxYc1sPrY5F3Jhaldgkm4EvgZ8MiJ+NH75noi4JOnNwJOSno+IpzYfO67YTgDcrDcm3ZJm1hV9yIlfq4SSlhlVXg9HxCPrr0fEpfG/l4GTjNJmmFnPjdLptJfQcFbqjEIKeBA4HxGf3/D6QWAw7hc7CHwY+NzMSlocT06yMpA0gbPQBiXnLA3JZ4t6pPuSvyFZ9ogs7KGwL8s4kWW3yLJRZBkdit9a9mczyfSQ/q5kmSWKIRZNf/maHafkZz1cq/5VvJaESqwNZ9tC2il9YPcAHwe+I+ns+LV/BzwPnBzVb+wBvhgRj82ikGY2X6NsFN1/hKwzCvlNyn9y3tVuccysC0ZTiXZABWZmu9EOaYGZ2e60UyLxzWyXWR+F7DpXYGZWyY+QCxZJqITaDjJOM04k4QtJyIbUcqaKJhksIM1ioSR8oRh+kS3OkWSVyKI5cqUDs5jq5GLZYUmrJfkYwFr1B3ItCaNYGc6uhbSeE7/rdnQFZmbNBLDqFpiZ9ZUfIc2sn+ovmbZQrsDMbIudlNDQzHYht8C6LJlgXVrmHUCFSeWlXPmjY+Y80bswehnJiGE2mkgyGjrIjivlsM8G+NLfmbb7ZLKCJEONWSHTHAHJ/S+MYg/XkgngMx2FdAVmZj0ViNWk8uwKV2BmVsl9YGbWT+FHSDPrKfeBmVmvuQIzs14KlI6AdoUrsO0qTYhOwxCaDddHg1CJ7HppHv3kfE1DLEqHzebXouWzNgilgTxUIp3MXbhe+mPOztcCd+KbWS+FO/HNrM+ygO6ucAVmZhU8mdvMeswtMDPrpQgYJgMSXeEKzMwqeRSyr7L89qXE7Fk4RDrsnqYvKO/LMlWUQj2yvPelzBEAq9sP2Rjt2v4vwDxDLJL0++m+tGWSrcOQ/DhL4Rdr2a2fcTaKPjxCTvy8SDok6euSzks6J+n+DfuOSvqepAuSPj3boprZ/Iw68etsi1TnD94q8KmI+HngvcDvSDosaQl4APgIcBi4T9Lh2RXVzOYpot5WR53GjqQ3SPqqpOfHDaZfnHTeiY+QEfEK8Mr466uSzgO3Aj8FXIiIF8YX/zLwUeC79b4lM+uyth4hNzR2PgRcBJ6R9GhEbK4rvgA8FhG/LmkvcGDSubfV5SDpduBO4GlGldhLG3ZfHL9WddxxSaclnV7h2nYuaWYLMBqFHNTagFvWf7/H2/FNp7uLcWMnIl4D1hs7r5N0M/BLwIOj68drEfHqpHLW7sSXdCPwNeCTEfEjVffQVjYoI+IEcALgZr2xZqPTzBap7uMhcCUijiT7qxo7d296z9uAHwD/TdK7gDPA/RHxk+zCtVpgkpYZVV4PR8QjGwpxaMPbbgMu1TmfmXVfhGptNdRp7OwB3gP8l4i4E/gJMHFgcGILbNzSehA4HxGf37DrGeDtku4AXgbuBT426XxzlS3ckaQKUfUq76PjCn+W0h9jKawBGoc25JkqCsclw+5ZdovsD7Gy8IuWNQ2xKIVEKF2cI/l8pOEQ5X1Z+EUpnGZtNTkm+ehMK6hdOdVRp7FzEbgYEU+P//9ValRgdT4T9wAfB94v6ex4OxYRq8AngMeB88BXIuJcjfOZWQ9Eza2G1xs74875e4FHr7tWxP8GXpL0jvFLH6DGgGCdUchvUmhgRMQp4NSkc5hZzwRES1OJImJV0npjZwl4aL2xI+kU8NsRcQn4XeDhcSX3AvCbk87tSHwzq9RmJH6psRMRxzZ8fRbIBgO2cAVmZpW2MQq5MK7AzGyLvsyFdAVmZlsF6boAXeEKrEoSfkEp/CIbP18qx2WkmSqyEItEqfM1/Tg2zCqRLgaSXa9l6XB66RanWSXKu7KFO/J9yeUKIS7D1ezzUT5fG/wIaWY9pdZGIWfJFZiZVXMLzMx6KdyJb2Z95haYmfWXW2Czlw2VNMjJ3rwYyWhRVsZsondW/jRpe/XwVDopOyvHanJccs7S9Zr+VLL7mI6YFSZKa9h0wnY2wTo5ZzJqWMp9nx8z4893cg+6ov8VmJm1z3FgZtZnjgMzs/5yBWZmveVHSDPrq2ycqCtcgZnZVqF8rmhHuAKrENny8IPS2HKSSD8LURiUh93z0Iws+Xp5V/FayT7tST4mDUIs0u8rufdEulhBcdegcFxWjux8ytZTGGaTr8vHDQq570uvAwyTe98Kt8DMrLdcgZlZb7kCM7NeciCrmfWZRyHNrL9cgZlZX7kF1mVZGIIaLGKfDP8nq9TnWSCyXPRJUVSKNsj6NJL8+22HWKRZJZJrpaENTUIshuVjlpJwiLWmoRLZOZerj1OScWKwMuM+qp3SBybpIeBXgcsR8c4Nr78IXGUUebQaEdtalNLMOiroxSNk3abGnwFHC/veFxHvduVltsNEzW2BarXAIuIpSbc3vYik48BxgP0caHoaM5ujLKljVzTo7LlOAE9IOjOupKrfFHEiIo5ExJFl9k15STObi53SAkvcExGXJL0ZeFLS8xHxVBsFM7PFUfRjFHKqFlhEXBr/exk4CdzVRqHMrANC9bYFatwCk3QQGETE1fHXHwY+11rJ+iYLy0iyF7CUfAAaLgbSaDGNpSQMoWmIRSnTRhLykC6AkpQxz9xRfVyWOSKSn9kg2RdJGEW2T6vVZRnsKf/USqEXrelBC6xuGMWXgF8GbpF0Efgs8HXgpEaxSnuAL0bEYzMqp5nNWR8eIeuOQt5X2PWuFstiZl0R/RiF3L2R+GaW2yktMDPbhVyBmVlf9aEPbMbDGGZms7OzW2DZkHyS6SELiSgNr2vQLB4mX+Ci2WIgpRCLNOQh2dcoO0dyvfRaaehIs4U2ivc4uYfZ+bIsFlpKPjtZiMVS9b4moRetafH0ko4CX2C0+s2fRsQfFd63BJwGXo6IX510XrfAzGyr8ShknW2ScaX0APAR4DBwn6TDhbffD5yvW0xXYGZWrf5cyFsknd6wbZ4XfRdwISJeiIjXgC8DH918OUm3Af8a+NO6RdzZj5Bm1ojYVif+lQnptG4FXtrw/4vA3RXv+xPg3wI31b2wW2BmVq29bBRV3Z7XHSlpPWHqme0U0S0wM9uq3WwUF4FDG/5/G3Bp03vuAX5N0jFgP3CzpD+PiN/ITuwKrCWRjZANmk30bpxLvzS61mACOCQ59mGUTLx40uqzxiAbxWs2gb3RBPFkFDKbzK1hUo7CaCKAVpPr7SmMQq4kH4LCMa1pbyrRM8DbJd0BvAzcC3xs4xsi4jPAZwAk/TLw+5MqL/AjpJkVrOcEm7RNEhGrwCeAxxmNMH4lIs4BSDol6a1Ny+gWmJlVazEOLCJOAacqXj9W8do3gG/UOa8rMDPbqgPpoutwBWZmlfowF9IVmJlVcwVmZn3lhIY7UWmid8MJz01z6TcKsWgwARwmhFhkzxkNFnyI7H4k9zjNpT8slCPJsZ9O5m4YfpElECiFWCgJy8gmek/NfWBm1ldiQtaQjnAFZmbV3AIzs77yKKSZ9ZcrMDPrJS+rZma95hZYhzXNl9/kUmmmima5+ZuEWKQZLDLJvYp0bYFCaEByP7LQBgZJORqEXzQKvYA8/GI1O277IRalXPmja802F0Mf+sBq3QFJD0m6LOm5Ta8flfQ9SRckfXo2RTSzhWgvoeHM1K3C/ww4uvGFbSbqN7OeaSudzizVqsAi4ingh5terpWoH0DS8fWE/ytcm6rAZjYHwSihYZ1tgaZ5iK5K1H9r1Rsj4kREHImII8vsm+KSZjYP64t6dL0FNk0n/sRE/WbWYz34bZ6mAquTqN/Meiodpe2IaSqwiYn6e6tJiEXDLArzDLGI7E9qcq1GC4hAMcNF+muR3HtlIRvZvSrc/0hTejTL3JGWcS0JESkc1zgsY1odGGGso24YxZeAvwTeIemipN/KEvWbWf/tmD6wiLiv8Hplon4z6z9PJTKz/urBI6QrMDPbqgOPh3W4AjOzaq7AzKyP1gNZu84V2DzMIMQio0GD3teGC4gwHJbLUQopyEIvkvM1yXwxKkfhPjb8uWQhG2kZh1k4SiGMIjtfFmLRgjRrSEe4AjOzrXoSB+YKzMwqOYzCzPrLLTAz6yt34ptZPwX5nOCOcAVmZpXcB7YTlf4qNV0IJBvKzzQIv2ic3aIc2ZCXo9SJ0jTjRKbJwiPZ95WEejTOipEpZcxIs2zMLozCcWBm1l8RfoQ0s/5yC8zM+ssVmJn1VR9aYLNd2tfM+imAYdTbapi0CLakQ5K+Lum8pHOS7q9zXrfA2tK0w3OOo5eRTNhuqvWJ4w3z9jeaVJ7JRjWz45r+PAujno1HZVvQVgtswyLYH2K0GNAzkh6NiO9ueNsq8KmIeFbSTcAZSU9ues8WboGZWbX1kchJG9yyvnD1eDu+6UwTF8GOiFci4tnx11cZrbNRuc7sRm6BmVmlbbTArkTEkWR/1SLYdxevK90O3Ak8PenCrsDMbKt20+nUXgRb0o3A14BPRsSPJp3YFZiZbSFANTvoa6i1CLakZUaV18MR8UidE7sPzMwqKaLWVsPri2BL2stoEexHr7vWaLTiQeB8RHy+bhldgZnZVrGNbdKpkkWwJZ2S9FbgHuDjwPslnR1vxyad24+Qi5b9BWt7CH2OE8dzSchDGirR4FJAZPntiwe1H8WZhkQUwkAiWz9gptqdC1laBDsi1iupS1T3laWmrsAkvQhcZfTxWp0wGmFmPdGHSPy2WmDvi4grLZ3LzLrA2ShGxoFtxwH2c2AelzSzaUSro5Az08YDdgBPSDpTEYE7ekPEiYg4EhFHltnXwiXNbOZa6sSfpTZaYPdExCVJbwaelPR8RDzVwnnNbIFqhkgs1NQtsIi4NP73MnCS0bwnM+u7+nMhF2aqFpikg8AgIq6Ov/4w8LlWSmbtfzjmnbe/0aVmETZQHaIwi5CNTKNwjmFy75vm368jgF2wqMfPACfH8S17gC9GxGNTl8rMFkrUjrJfqKkqsIh4AXhXS2Uxsy5Z634TzJH4ZrbVLnmENLMdasc/QprZDuYKzMz6afEhEnW4AttN5v2BbLSYxhxDNrJQiSYhDxO1G5uhmHEYRQ+mErkCM7NK7gMzs/5yBWZmvRRAo2SV8+UKzMwquBPfzPrMFZiZ9VKQTyTvCFdgNjs9+Ate1v1f3jQMZPqzzzWkpSlXYGZWrQd/gFyBmdlWHoU0s15zC8zMessVmJn1UkRxtfAucQVmVqUHrY+Z68E9cAVmZtVcgZlZP4VHIc2spwLCgaxm1lueSmRmvRThZdXMrMfciW9mfRU9aIFNvXKBpKOSvifpgqRPt1EoM1u0cULDOtsCTVWBSVoCHgA+AhwG7pN0uI2CmdkCrU/mrrMt0LSPkHcBFyLiBQBJXwY+Cnx345skHQeOA+znwJSXNLNZCyB6MJVo2kfIW4GXNvz/4vi160TEiYg4EhFHltk35SXNbOZinNCwzrZA01ZgVStrdn/owswmirWotdVRp6+8SX/6tBXYReDQhv/fBlya8pxm1gUttcDq9JU37U9XTDGKIGkP8DfAB4CXgWeAj0XEueSYHwB/t+GlW4ArjQvRHpfjel0oRxfKAP0sx7+IiDc1vZCkx8bXq2M/8M8b/n8iIk5sONcvAv8xIn5l/P/PAETEH27nPVWm6sSPiFVJnwAeB5aAh7LKa3zMdTdV0umIODJNOdrgcnSvHF0ow24tR0QcbfF0VX3ldzd4zxZTB7JGxCng1LTnMbMdq05feaP+9KkDWc3MJqjTV96oP70LFdiJyW+ZC5fjel0oRxfKAC7HtJ4B3i7pDkl7gXuBRxu8Z4upOvHNzOqQdAz4E/5/X/l/Gr9+CvjtiLhUek96XldgZtZXXXiENDNrxBWYmfWWKzAz662FVWBdySMm6UVJ35F0VtLpOV/7IUmXJT236fW53pukHHO5N5IOSfq6pPOSzkm6f8O+ed+LrCzzuh/7Jf2VpG+Ny/AHG/Z14vemMyJi7hujUYbvA28D9gLfAg4vqCwvArcs6Nq/BLwHeG6R96aqHPO8N8BbgPeMv76J0fS0wwu6F5VlmfP9EHDj+Otl4GngvV36venKtqgW2Ot5xCLiNWA9j9iuEhFPAT/c9PLc702hHHMTEa9ExLPjr68C5xlNLVnEvSiVZW5i5Mfj/y6Pt8C/N1ssqgKrlUdsTgJ4QtKZceLFRdvV90bS7cCdjFodC70Xm8oCc7wfkpYknQUuA09GxMLvRxctalGPLuURuydGQXRvBp6U9Py4RbIou/beSLoR+BrwyYj4kaSF3YvNZRm/PLf7ERFD4N2S3gCclPROuvXZ6IRFtcA6k0csIi6N/70MnGTUTF+kXXlvJC0zqjAejohHxi8v5F4UyrKQz0pEvAp8AzhKhz4bXbGoCqzRvKe2SToo6ab1r4EPA8/lR83crrs345bWg8D5iPj8hl1zvxelssz5frxp3PJC0g3AB4Hn6chno1MWNXoAHGM0wvN94N8vqAxvYzSS8y3g3LzLAXwJeAVYYfTX9bcWcW+qyjHPewP8K0aPQt8Gzo63Ywu6F5VlmfP9+AXgr8dleA74Dxv2Lfz3pkub50KaWW85Et/MessVmJn1liswM+stV2Bm1luuwMyst1yBmVlvuQIzs976f8JUFNNXZSlzAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "im = ax.imshow(np.transpose(v.v()), origin=\"lower\")\n", + "fig.colorbar(im, ax=ax)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "14d6cd34-62a6-4c75-a39b-5aedbbc97c93", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "e = v - true(mg.x2d, mg.y2d)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "13e70155-7822-48cc-90ac-08199b0634c2", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.000617617264475653\n" + ] + } + ], + "source": [ + "print(np.abs(e).max())" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "aee1546b-2652-4b65-a02e-c3b507973676", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVIAAAD6CAYAAADz7c/YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfy0lEQVR4nO3dXaxc1XUH8P9/5l7b+NoNCeYrtluIhNJaSEBkmSBXVUgJdVAkXgGJp0R+KRJIeUHhIUqkSu0Lah94iFUsoooPoYJVlN5ijESEkCKwQSZgjCvHJeJiKseA8ef9mJnVhxnTuXf2WrPPmbln5tr/nzTCzD77nD3nzl33zNlr1qaZQUREyquNegAiIiudAqmIyIAUSEVEBqRAKiIyIAVSEZEBTVR9wA3fqNsNmycL9zOkswuinAOvT79+rSCTwWtpgv7+gram+X/LmsHfuYbVC++zER0raAv7taJ9pl93K+hjLf9coeU3IejHpvN8sL+wzdkfANSa/nuHjaCt6Ryw4R/MGg237Qy+OGlmV7sb9PF3d07ZZ58HL7TL27+f22tmO8oe61JQeSC9YfMk3tq7uXC/pqXfaA34P+wF89tmwzb/DT/nNJ1p+X8czpnfdqq11m0707zCbfusuc5t+7KR3udnC1Nun8+Dti/m/DGemvPHeHp2dfL5cxfSzwPA/Pngj+wF/4/HxNmg7Vw6yE6c8w+16oz/HojaVn/pv69WfTHvtk18cT7dcPKU26f5pz+5ba/av//Rbcxw8vMm3ty7KWvbyev/sGGQY10KKg+kIrISmHvxIr0USEWkhwFohTfApJsCqYgktcKb0tJNgVREehgMC/pon02BVER6GICmPtpnqzyQLqCFE830dGkzmC335kKjv5kLwftgNkgfitvSp+xcy5+JPm9+2+nmGrfty6Y/k/5lMKP/ZSPddtp5HgDONVa5bRca/kz6XNM/VwtOW6sZpC87KVMAwKAt5HULhtGaCFLWgsSC5uogjWyt/+tWW0i/D2qN9W4f/8wDOBE15tE90nx9E/JJriH5Fsl3SR4i+Yuuth0kj5A8SvLR5R2qiFTF0L6wyXlI3hXpHIDvm9lZkpMA3iD5XwD2A3gCwA8AzADYT/IlM/tg+YYrIlXRHdJ8fQOptQuWnu3872TnYQC2AThqZscAgORzAO4FoEAqssIZTPdIC8j6rj3JOsmDaN952WdmbwLYCODjrs1mOs+l+u8keYDkgc8+0985kXFn1p5jyHlIZiA1s6aZ3QpgE4BtJG9G+hZ+8rSa2S4z22pmW6+6SnVSRMYf0cx8SMHqT2Z2CsBvAexA+wq0+0vzmwAcH9bARGR0DEDL8h6ScY+U5NUAFszsFMkrANwF4J/Qnmy6ieSNAD4BcB+AB/rtb8GA/w1SZjwtJ/Ul+ou4EKYx+Tks80E/L5VpNipaEqRGnWn56U9no9SoKJWp6RQLCVKczgdtc03/bdIIfpYtpyKTRWlMJX8xw13W0jtt1YOKXcFvRvCjRmNNVIXKP1dspc9/VCetNhH8Hg0h/UlXm/lyZu2vB/BrknW0r2CfN7PfAADJhwDsRTulbbeZHVq2kYpIZdoJ+QqkuXJm7X8P4DanbRrA9LAHJSKjZQAWgjq0spi+IioiPQwMC4vLYgqkIpLkzUtILwVSEemhe6TFVB5Im1bD58300hVlfnALThGRfvubdWZJ2/ssPtsfFi0JjnXemWHv18+bmQeAs84M/Gww++4VGAGAhWCNpUZQgCRam2nonJl5ADBndj74MYcz883Vwcx8uA5U8fNhNf89MLGqePZLPobreMliuiIVkR7tCvkKpLkUSEWkhxnDfGpZTIFURJKiZcRlMQVSEenRnmzSR/tcCqQikqDJpiIUSEWkhyabiqk8kDZQx2fNdYX7ealMreCvZnSzPEqbKpP+NBfky8TpT8FaSWGbf7xZp20+qMQx3/JfczNI24mStr1VKMquTuEVH2m3BQVInJfGIMWJ3iJhAJpRSd2SSexWS5/jaO2oVrA+1DA0lZCfTVekItLDwPBiQxbTmRKRHppsKkaBVER6GKiP9gUokIpIkiab8imQikgPMyj9qQAFUhHp0Z5s0ldEc42k+tMpp/pT2M/5mBGl30Q3y6M3SdTmpTlFazbNBWlHZVOconWU5p22qMJT6RQnt8XH6NZbyTarF6/+1Jr0+9BZb2ogQYqW95ZrTgbv76AK1TBosimfrkhFpIeBKuxcgAKpiCTpijSfAqmI9Giva69AmkuBVEQSqKVGClAgFZEe7eWYNWufS4FURHqYUR/tC6g+/QnEmdaawv28H2r08SN6I5RNf1pwSglFKU5RZagwxalktaaG87q954G40o+VbHMxSDsKKzwFyVZOihMAmJPmFI29FZaoCl5zmKLlt3lVnmr+GodoLvhtw6CE/Hx9zxTJzSRfI3mY5CGSD3e1fUTyPZIHSR5Y3qGKSFXa9UiZ9RgGkjtIHiF5lOSjRbfx2kr2KRzXcq5IGwB+ambvkFwP4G2S+8zsg077nWZ2MudgIrJSVFchn2QdwBMAfgBgBsB+ki91xZhwG68NwJGifcrGtb5nysw+NbN3Ov8+A+AwgI25BxCRlaed/sSsB4ANJA90PXYWPNw2AEfN7JiZzQN4DsC9Bbbx2sr0KaXQPVKSNwC4DcCbnacMwCskDcCvzGyX028ngJ0A8PXri98fFZFqFfyu/Ukz2zrA4TYC+Ljr/2cA3F5gG6+tTB8gM651yw6kJNcBeAHAI2Z2uvP0djM7TvIaAPtIfmhmry/t2xnILgDYfPOflVxoQkSqNMwyeiRfBXBdoukxpKfolsaJaBuvrUwfIDOudcsKpCQn0Q6iT5vZi18d1ex4578nSO5B+3I5PKCIjL92Gb3hJeSb2V1eG8k7AGzuemoTgONLNpsJtvHayvQpFdf6BlKSBPAkgMNm9njX81MAamZ2pvPvuwH8st/+WlbD+WaQ0+H1c2YHoxvi0Yyil8YElEuNihaWi48VLN4XpTiFbcUrZcVtblPIq/IU/nqGKU5BalQ0SKcpfl1BOlWUolUixQkAak4WXLQIH5vL+82jCouW7AdwE8kbAXwC4D4ADxTYxms7UrRP2biWc0W6HcCDAN4jebDz3M8AfAhgTzvOYgLAM2b2csb+RGTMtas/VTNrb2YNkg8B2AugDmC3mR0CAJLTAH7S+aid3KZP/0J9SH4LJeJa30BqZm/A//N8S7/+IrLytL8iWl1CvplNA5hOPH9Pv2369C/Ux8yOoURc01dERSRBXxEtQoFURJKG9a2ly4ECqYj0GPas/aVOgVREkvTRPl/lgbQF4nzLX/DN7ef8dYyWQyib7hOlPzWctijFaS5oi9KmvDSm9jiiSk7ptlKVmvpgVMnJaWOt5fcJFogL05/M36d575HgUM0ghtQafpsF42ew2F7LSXOi/7KWZ4G+Dq3ZVIyuSEWkhyH+Yy2LKZCKSJI+2udTIBWRXqaP9kUokIpIj4uFnSWPAqmIJOmKNF/lgdSM4VpEnjI/1Gj2PV7PKZgtd2bgy/Rpt0VZB+UyEobZp59oj96sfa3kukytejCFHWRvGNL9zKuqAgDB7HsrWB8Kwcx8VGTEnZ2PCqssY0HKi4WdJY+uSEWkh4HhH3lZTIFURJJ0jzSfAqmI9DJ9tC9CgVREeugeaTEKpCKSpECaT4FURHoYiKYmm7JVn/4EPy2pXEpPuTWbohJhUbqSt89oHNF3lqPXXLbfsK8kojoitaBoiZfmZBateVQuxakVVvdIvwCL0piCVKWwX1BIxKJ8Ja9phGvuarIpn65IRaSHabKpEAVSEUlajrKLlyoFUhFJUNGSIhRIRSRJV6T5FEhFpIcZ0FzGCvyXGgVSEUnSrH2+EaQ/0U0vKvODi9KYylZPisbhrc1UPlWpXK7esD92RWlMZdZlitrC6k9BahScKk7tYwXn2MnfsihVKUpxCisyBfsskcoUFahaTgZ9tC+i728xyc0kXyN5mOQhkg93te0geYTkUZKPLu9QRaQ67cmmnIfkXZE2APzUzN4huR7A2yT3ATgC4AkAPwAwA2A/yZfM7IPlG66IVKXMFfTlqm8gNbNPAXza+fcZkocBbATwNQBHzewYAJB8DsC9ABRIRS4B+mifr9ANOpI3ALgNwJtoB9OPu5pnOs+l+u0keYDkgQtfzJYcqohUpT1rX8t6SIFASnIdgBcAPGJmp5FeZSL5YcDMdpnZVjPbesXX15QbqYhUyizvIZmz9iQn0Q6iT5vZi52nZwBs7tpsE4Djwx2eiIyKPtrn6xtI2c4reRLAYTN7vKtpP4CbSN4I4BMA9wF4oN/+2tWfin8cKJMmNOwUp2ifZasxReMYfhUn//IhSiMrmxpV99pqUYUnX5TiZEGJKjppTuFbquS5j4JPdK78TsWraw2DgQqkBeREp+0AHgTwfZIHO497zKwB4CEAewEcBvC8mR1axrGKSIUs8yF5s/ZvwFl118ymAUwPe1AiMmIWf2FBFtNXREUkSR/t8ymQikiSZuTzKZCKSA99174YBVIR6WUonblwOaq++pNx6KlMbp+Si9+V3edKFr2qMMUpTMFx0pyCb8MwSI0K05+i9K0x+fJNVMmp5rzu6OcyUW8ONqA+9NE+35i8xURkvBDWynsM5WgZleSibby2Pn12kzxB8v2iY1lKgVRE0ipKJCVZR7uS3A8BbAFwP8ktudt4bRn7fQrAjqJjSVEgFZFe1r5dkvMAsOFiUaLOY2fBo21Dp5Kcmc0DuFhJLncbry3cr5m9DuDzEmPpockmEUnLv9o8aWZbBzhSqpLc7QW28dpy9ltmLD0USEXEMbyJVZKvArgu0fSYc6ClYTzaxmvLrlCXeRzXSALpMItxlJ2ZL7ue0zD7DCKaSa85P/dW1Cdoi5coClrdUxzMzEez72UzLUr0W441rIK6Ku4MvFv4BcDkMs/aBz+mwszsLq+N5B3oX0kuqjbntZWpUFeqqp3ukYpIr4t5pDmPwX1VSY7kKrQryb1UYBuvLWe/ZcbSQ4FURJKqKuwcVZIjOU3ym9E2Xlu/CnUknwXwOwDfJjlD8sdlq9rpHqmIpFWYkO9VkjOze/pt06d/1Of+IvuKKJCKSJq+IppNgVREksoU9L9cKZCKSC8joMLO2RRIhyRKlymbGjXsfYYpTiWPVY+KjLj9wsWS/JYhf9SMU5WGvE4VgIngXHmpTJM1P8VpTX3BbRsKXZFmUyAVkTQF0mwKpCKSpkCaTYFURHqpsHMhCqQikqRZ+3wKpCKSpkCaTYFURJJ0RZrvkg6kUSpKXNpmuJWhIl6lJiCubBWmMjlfgI72VzYVKDofbr8gDSgSpT9F4/dEryt670QpX2GKU4lUpjX1htvnimVPf9I90lxZgZTkbgA/AnDCzG7uev4jAGcANAE0BizuKiLjYkjLiFwucqs/PYUla5t0udPMblUQFbnEVLRm06Ug64rUzF4neUPZg3TWcNkJAFPXTZXdjYhUiEMs7HypG7QeqQF4heTb0YJXZrbLzLaa2dY1V64Z8JAiUgldkWYbdLJpu5kdJ3kNgH0kP+yszCciKxhNs/ZFDHRFambHO/89AWAP2kuZisiloLqlRla80lekJKcA1MzsTOffdwP4ZU7fKO3E46XZlE0fKss9XnCo+PX6N6JqwU6jtCMvFSg6V5FwgbuAN8YotSg6UpkUJ8A//9HPJVpYrh7cPFwV9IuqNXlt6ybm3T5T9Tm3bSh0RZotN/3pWQDfA7CB5AyAnwN4DcAekhf384yZvbxM4xSRiumjfb7cWfvk2iYAbhniWERkXJhm7Yu4pL/ZJCID0BVpNgVSEUlTIM2mQCoiSbpHmm/QhHwRkcte5VekpKFW6i52iZgf5tJEHYtXhipVBalfv+AFlEkh86pCAQjPRyuqDBUczxvjsCtoRccC/LSpspWaVkVtQbWmqTCVKd02NeGnOH2tfsFtGwpdkWbTR3sR6aVZ+0IUSEUkTVek2RRIRaQHocmmIhRIRSRNgTSbAqmI9FL1p0JGEkjjtZS8PumZ0mY4AxzMKQ95Rj86kY1gHNHMcSO42R8VNPGEs+/B+QizDoLjeWsslck46CcqaOKd46j4SPRziWbmo8Ik0RpL3ux8NDP/jYmzbttQaLIpm65IRSRJV6T5FEhFJE2BNJsCqYj00jIihSiQikiSPtrnUyAVkTQF0mwKpCKSpK+I5qu+aAmACSeVqYxakFoUpftEaVMt8/fppR35iS1ALUrRKrnWU61EAZJof9H5KJuu5P0eemlR/UQpTuG5ci6tyvTp128yWs+p5qdNra2li5asr8+6fa6sn3fbBqZ7pIXoilREehB90qllEQVSEUnTFWk2FXYWkSRa3mMoxyJ3kDxC8ijJR4tu47X16bOb5AmS7y95/iOS75E8SPJAzvgVSEUkzTIfAyJZB/AEgB8C2ALgfpJbcrfx2jL2+xSAHc6w7jSzW81sa85rUCAVkV6dws45DwAbSB7oeuwseLRtAI6a2TEzmwfwHIB7C2zjtYX7NbPXAXxecKxJukcqImn5V5snc6/cHBsBfNz1/zMAbi+wjdeWs98UA/AK22kivzKzXf06jCD9yTA5xPSnZun1kMqlTS206snnJ6NSOcF1f8PZH9AnPSdK8nPSt8IqTtHNrmVYY2ncRec3So2aCNZzWh2kP62ppRPo1tb8NZvWL/OaTcP8ZhPJVwFcl2h6DOkEgaVHj7bx2nL2m7LdzI6TvAbAPpIfdq5eXVmBlORuAD8CcMLMbu56fgeAfwFQB/CvZvaPOfsTkRVgiIHUzO7y2kjeAWBz11ObABxfstlMsI3XFvWJxnq8898TJPegfYsgDKS590ifwpKbsjk3iEVk5apw1n4/gJtI3khyFYD7ALxUYBuvLWe/i18zOUVy/cV/A7gbwPtRHyAzkDo3ZXNuEF8c3M6LN6IvfOF/VBGRMWFofzUt5zHoocwaAB4CsBfAYQDPm9khACA5TfKb0TZeW9Sns+9nAfwOwLdJzpD8MYBrAbxB8l0AbwH4TzN7ud9rGOQeafaN3M7N2l0AcO2WbyjNV2TMVb34nZlNA5hOPH9Pv2369I/63O8M55aMIS8ySCAteyNXRFYC/TZnGySQlrqRKyIrA6PCOLLIIIH0qxu5AD5B+0buA/061WhhGoin6dzOnQhu0kTpT1HbggUpSc4qcQtBxajwPlJ0l7pkP2/RvGjBvCg1yoJfqFaFpS2iqlHRqWox3Y8lK4CVVQ9G6aUETkXpTzW/MtTAVP2pkKx3S+qmbL8buSKyslX5XfuVLuuK1LspG93IFZGVTYWd8+kroiKSpqvNbAqkItJLH9sLUSAVkTQF0mwKpCLSo+qE/JVuJNWfyqQ/eaL0m2aQwhKmOAXpPm6/ilOcWq3iqUzhgnnLwEtXKpuWFim1QF+UQha8d6K2smlTXrWpqFLaFNML5g0LW4qkuXRFKiK9lEdaiAKpiCQp/SmfAqmIpOmKNJsCqYgkabIpnwKpiPQyACpakk2BVESSdI80X+WBtEbD2nrxKvleWkmU4uRV/QGASfPTSsqkRkWVfWrB/kLBG7kVvMu91x0u6Eb/PA5vqcK2KMUpugaK+jVLtEXVpKIF7ubp/9rM1f2fdfS+8tSDn9maIS4iuZTySIvRFamI9DLTR/sCFEhFJElXpPkUSEUkTYE0mwKpiCTpijSfAqmI9DIATUXSXNXP2sOwtla82II7ax8ULYkKSET9Flr+afFmXueCU7kcxUKiYi3e7HaUxVB6HCUKkIQz7C3/Zxau2RSc4prTrVnzZ8QbwTgiEzV/Jv3cxGq3ba41WfhYyzlrD+iKtAhdkYpImmbtsymQikiSrkjzKZCKSC+V0StEgVREehAANdmUTYFURJKoe6TZFEhFpJc+2hdSeSCto4X1tQuF+7WcBXaioiXNYFGeqIDEQpBW4vWLikvMlkht6SdM7XJScKK1hkqtedSHl+YUpTg1msHrCvqVWeuJwWueCFKjovStKC1tVZAadeXEFcnno/dOuQStXPqufREDB1KSHwE4g3aRoIaZbR10nyIyepq1zzesK9I7zezkkPYlIuNAV6TZKvloT3IngJ0AcNU3V1VxSBEZhGnWvohh3GYxAK+QfLsTMHs3MNtlZlvNbOv6rw//fqGILAPLfMhQrki3m9lxktcA2EfyQzN7fQj7FZERUvpTvoGvSM3seOe/JwDsAbBt0H2KyBi4WCW/30MGuyIlOQWgZmZnOv++G8Avoz51tnBl/XyyzUtxikTpT1GK03yU/mTFqz9NBilT0fo/UdpUJEp/8lKBGtGaTcEYozShuCKTl/7k94mqLjUa/s+sFezTbQvGzpr/mucW/PfHfDDGyNREuiLaqdVr3T7LujadLfcBLi2DfrS/FsAetsuzTQB4xsxeHnhUIjJShOmjfQEDBVIzOwbgliGNRUTGSUuXpLmW98sRIrIyXfxon/MYApI7SB4heZTko0W38dqC5zeTfI3kYZKHSD5cZCxLKZCKSBLNsh4DH4esA3gCwA8BbAFwP8ktudt4bX322wDwUzP7KwDfBfD3GX1cCqQikpY/a7+B5IGuRzKfPLANwFEzO2Zm8wCeA3BvgW28NrePmX1qZu90/n0GwGEAGzPH0kPVn0QkoVBq08kBa2xsBPBx1//PALi9wDZeW85+QfIGALcBeBPtzKO+fZaqvvpTkP4UKbP4XZTGVDb9adbS38yabPnpT/Uo/ak1/K/MehWIFoKUqWjRNjTLvU28ak2Npn/uoxSnxoLf1gyqRlkj3WZBylRkIUqNmvR/nvONYIFE5z1y7erTbp/ZteVSrbIMeRVRkq8CuC7R9BiQfMMuPXi0jdfWd78k1wF4AcAjZnaaTK4Q2fdE6IpURJKGmf5kZne5xyHvALC566lNAI4v2Wwm2MZri/qA5CTaQfRpM3sx4zgu3SMVkbTqvtm0H8BNJG8kuQrAfQBeKrCN1+b26Vx5PgngsJk9XnAsPXRFKiK9DECrmoR8M2uQfAjAXgB1ALvN7BAAkJwG8JNOPY/kNn36J58HsB3AgwDeI3mw89zPzGw66ONSIBWRhGq/R29m0wCmE8/f02+bPv29599A+h5qeByPAqmIpOkrotkUSEWklwFo6iuiuSoPpBNo4aoSi995aU7RomfzJRe/81KcorZZ+mkvk2y4bbWoIlPJylDea1tdC6phBYu9lV0Yz6sMFf3MGk6qEhCnP9l8MG/q7JMLQfWnZpAaFZyOoHgVzq7231f/M5tuWz3hv3fuXL/OP9jADDAF0ly6IhWRNH20z6ZAKiK9Kpy1vxQokIpImq5IsymQikiaAmk2BVIR6WUGNIP6C7JI9bP2BK6uF58NbDp/HaMf9ULwB3U2KOAxZwtBv/QpO0e/TzRrXw9mxOtB1dxozaYFZ3a+EczaX6A/ozwRzOiXERUYaQUFTWwhmBIPZu3rF9JttXl/Zr4+57fV0ssrtYUXcUFGwlT6db93YZPb5/Ur/zI41pFoIHl0RZpNV6QikqZAmk2BVEQSTLP2BSiQikgvA0wJ+dkUSEUkTV8RzaZAKiK9zLQccwEKpCKSpsmmbJUH0knUcE19qnC/pnO/phEkQC1Y1Ob/tZ0N0p/OO21rgz5rgtSoVYzWeip3RdB0UqOiQi0XmkH6UzTGEqlRXjETAGhFxUKC9KfabJT+lN7nxDn/WBPBsmKT5/wAE7VNzPptNSdXb2Gd/yv6b2f/2m0D/iNoy2O6Is028FIjJHeQPELyKMlHhzEoERm1zGVGdNUKYMBASrIO4AkAPwSwBcD9JLcMY2AiMkIXi5bkPGTgj/bbABw1s2MAQPI5APcC+KB7I5I7AewEgD/fqNuyIuPOAJi+Ippt0I/2GwF83PX/M53nFjGzXWa21cy2Xn3VMq7FLSLDYZ3CzjkPGfiKNHW3Xtf6IpcA08f2bIMG0hkAm7v+fxOA4wPuU0TGga42s9EGmHUjOQHgvwH8LYBPAOwH8EC0DjTJPwH4Y9dTGwCcLD2I4dE4FhuHcYzDGICVOY6/MLOryx6I5Mud4+U4aWY7yh7rUjBQIAUAkvcA+GcAdQC7zewfCvY/YGZbBxrEEGgc4zeOcRiDxiE5Bp5CN7NpANNDGIuIyIo0cEK+iMjlbhwC6a5RD6BD41hsHMYxDmMANA7pY+B7pCIil7txuCIVEVnRFEhFRAakQCoiMiAFUhGRAY0skI5LHVOSH5F8j+RBkgcqPvZukidIvr/k+UrPTTCOSs4Nyc0kXyN5mOQhkg93tVV9LqKxVHU+1pB8i+S7nTH8oqttLH5vZAkzq/yB9reg/gDgWwBWAXgXwJYRjeUjABtGdOy/AfAdAO+P8tykxlHluQFwPYDvdP69Hu2vHW8Z0blIjqXi80EA6zr/ngTwJoDvjtPvjR6LH6O6Iv2qjqmZzQO4WMf0smJmrwP4fMnTlZ8bZxyVMbNPzeydzr/PADiMdjnGUZwLbyyVsbaznf+d7DwM+r0ZW6MKpFl1TCtiAF4h+XanAPWoXdbnhuQNAG5D+ypspOdiyViACs8HyTrJgwBOANhnZiM/H+IbVbn6capjut3MjpO8BsA+kh92rtBG5bI9NyTXAXgBwCNmdprkyM7F0rF0nq7sfJhZE8CtJK8EsIfkzRiv94Z0GdUV6djUMTWz453/ngCwB+2PT6N0WZ4bkpNoB66nzezFztMjORfOWEbyXjGzUwB+C2AHxui9IYuNKpDuB3ATyRtJrgJwH4CXqh4EySmS6y/+G8DdAN6Pey27y+7cdK48nwRw2Mwe72qq/Fx4Y6n4fFzduRIFySsA3AXgQ4zJe0N6jeSjvZk1SD4EYC/+v46pWwx6GV2L9scmoH0unjGzl6s6OMlnAXwPwAaSMwB+bmZPVn1uUuMA8BqqOzfbATwI4L3OfUEA+JmZTY/gfZIcC9qBrKrzcT2AX7O9Sm8NwPNm9hsAGJPfG1lCRUtERAakbzaJiAxIgVREZEAKpCIiA1IgFREZkAKpiMiAFEhFRAakQCoiMqD/A8PLTAMvjUPOAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "im = ax.imshow(np.transpose(e.v()), origin=\"lower\")\n", + "fig.colorbar(im, ax=ax)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "cc53046f-d3f2-46b6-ba1c-9f80eca9b10c", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8.044915834004946e-05" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "e.norm()" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "923f23f0-8879-4a45-a6d7-0443a11c393b", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.9832702148195685e-12\n" + ] + } + ], + "source": [ + "res = mg._compute_residual(mg.nlevels-1)\n", + "print(mg.grids[mg.nlevels-1].get_var(\"r\").norm())" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "cd76edcc-228f-4448-9fb4-ad7d434a6c13", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.1026179561330153e-13\n", + "9.578225232788658e-14\n", + "7.273262526067396e-14\n", + "6.897528548269955e-14\n", + "6.897184446661625e-14\n", + "6.756122004635183e-14\n", + "6.733362509180792e-14\n", + "6.86361322648682e-14\n", + "6.754540465868575e-14\n", + "6.733362509134965e-14\n" + ] + } + ], + "source": [ + "for i in range(10):\n", + " mg.smooth(mg.nlevels-1, 100)\n", + " mg.grids[mg.nlevels-1].fill_BC(\"v\")\n", + " mg._compute_residual(mg.nlevels-1)\n", + " print(mg.grids[mg.nlevels-1].get_var(\"r\").norm())" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "cbf53be4-e327-4a26-bf71-27e954449306", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "r = mg.grids[mg.nlevels-1].get_var(\"r\")" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "e608eabf-2253-4bfa-8b64-6fde3f1f60af", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEJCAYAAADFK6HzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeYElEQVR4nO3de6wc130f8O939z55L58iKQkkJUqJ7Ih2a8tl1MRyU8sPlVYMqyhQQ1ITBIELIkXVymjQRK7RGjFgGAUKwwWqtiEqwg38EFzbRASb1SO1FFWoK5GyZUUUqZRiqOjyUrkiKZnP+9r99Y9d2le88/vdubOz8yC/H2BB7pw7M2dnd8+emfOb36GZQUSkbhplV0BEJAs1XiJSS2q8RKSW1HiJSC2p8RKRWlLjJSK1pMZLRGpJjZeI1JIaLxFJRPJGkg+R/E73+c0k/yvJ75D8Z2XXT42XyGWG5G6SUyRfSijbQfIVkodJPhBtx8yOmNlnFjw/aGa/B+DTALbnX/PlUeMlcvn5GoAdly4k2QTwIIBPANgG4B6S20j+LZLfv+SxMWnDJD8F4BkA/6t/1U9noOwKiEi+zOxpklsTim4FcNjMjgAAyYcB3GVmXwbwyZTbfgTAIyR/AOCbOVU5EzVeIjVE8nYANLMfdp9/BsDjZvZ6sNomAAvLJwD83WAfVwH4EoBbSH4OwI8A/CMAwwD29vYKeqfGS6SGzOxJkp8m+SEAWwH8ZImGCwCYtKlgHycB/N4li59aTj37Sde8RGrKzL4N4B8DWGFmP06xygSALQuebwYw2Y+6FUGNl0hNkfw0gP8B4DzJD6RYZR+Am0jeQHIIwN0AHulnHftJjZdIDXWveZ00s2fM7OvoXJe6rlv2LXSuT72b5ET3ehjMbB7AfQAeA3AQwLfN7EA5r6B3VCZVEakj9bxEpJbUeIlILRUeKtEcH7OBdeuK3m3/tYOypAHqixrBaXuwHpv+ejaX/JvElr+9wfE5f3tTg25Zc+OsWzZ9YSi5YCA4WOa/aM74ZU2/GjDvJzo49K2VQR3bfj0a035Ze9jfYWMweX+rh6fddc4dW+GWnf3ZsRNmtsH9gyX8g9vH7OSp4AOzwPMvzjxmZosi+vut8MZrYN06XPuH9ycXRpffogbAW6Xlr2RRoxFt0/ngRh9aG/D31RoPviRBAzW4csYtm39zNHH50Kmmu861Hzzmb++/XOOWrfwXfmjRoRevS1ze3OB/IVuzfh1HXh326/Gaf6zmnO84g0P/9m/4dWxf8L82q172G/rT7553y1ZsOJe4/JM3+tfTn/38r7pl//sHf/iaW5jCiVMtPPvY5lR/O3jtq+t72VdWClIVkQSGlkWnE+VT4yUiixiAdngqVD41XiKSqB1eyC2fGi8RWcRgmNNpo4jUjQFo6bQxQZ7HJNhWMOrujhoutZ67jj9IFg6RD77lr9geDMIhTo65ZRxJXi/a3tFXr3bLGn8vCFH40fVu2egZZ2R2Ihjif5cfsjG/wq//1K/7w/rDU8kf85mr/dG/X3rIf81Tt/gjim5YBoCRN/yv2y/dfDJx+XcO3OKu0/6H/r7wg6Aspapf81oySJXkCMnnSP6U5AGSf7SgLHVKWRGpDwPQMkv1KEuantcMgI+Y2VmSgwCeIfk/0blD/UEAH0cn1cY+ko+Y2cv9q66IFKXaV7xSNF7WuXP7bPfpYPdhcFLKAlDjJVJzBqv8Na9U9zaSbJJ8AcAUgCfM7Fkkp5Td5Ky/k+R+kvtbZ88m/YmIVIgZMJfyUZZUjZeZtczs/ehkXryV5HuxjJSyZrbLzLab2fbm+HjmyopIUYhWykdZlpVVwszeRieH9Q5cZillReQXDEDb0j3KsuQ1L5IbAMyZ2dskRwF8DMC/x4KUsgCOoZNS9t5Uey2qsWYQapCxEubcLO0tBwDO+fuaW+sP8TfP+b8tDT+iAPPONgfO+W/38N/4ZXNBhgULPkFzq5OPyfwKf3tDQTjB+PuTwwkA4K23/B79/Hjycbz2Sf/4tgf9spV3vOGWnXzODzm59v/4oRmv/ezGxOWv/uv/7K7zy0/+rluWhzJ7VWmkGW28FsB/705Y2UAndez3AYDkxZSyTQC765xSVkR+oROkWvPGy8xeBJAYKWdme1GB+dtEJF8GYC6KuK0A3R4kIosYiFbOiZa7Z2/7ARwzs1QzdEfUeIlIonaW++Ri96Mza9GqPDZW7X6hiJTi4jWvlKES6y/GcXYfOy/dHsnNAH4TwH/Lq47173ll/XUIUjNjfvnbjFJON2aDG5un/Ruzo0sOs5v84cahyeQbh+dW+6N80YioBTd0z68KUhu/mpzDvjGT7Wb02T/3sw0PrAvSQK9LHn2dWePXoxmko57e448oDo/6x3H4zfNu2bmPJndG3vOjf+Kus2XjKbfsr9yStIhW+mteJ8xs+xJ/81UAfwBgZS+1Wkg9LxFZpJNJtZHqsRSSnwQwZWbP51nH+ve8RCR3ZsRslOdpeW4D8CmSdwIYAbCK5NfN7Ld62ah6XiKSqA2meizFzD5nZpvNbCs6wew/7LXhAtTzEpEEnQv21e7bqPESkQTLumCfmpk9hc790T1T4yUii1y8YF9l1Wq8cp4xO+u+opmU3RuRgxulW6PBTdvRjOrBZ4dn/Yups1c5N2af8dcZPO0f4Nlt/hB/80jy7NwAMLfSuYk9eF3RNeKZv5M8qzQAtIO8+MNTyRsdOuO/L5N/3z8eG9895ZadC27MPv6h1W7Z9Y8mz4D+V7/jzxJ+9C3/2OehlX+Qaq6q1XiJSCUYiLkoZUgFVLt2IlIKXbAXkVoyUKeNIlJPumAvIrVjhr6ESuRJjZeILNK5YJ/b7UF9UZ/GyxvVDvLUh4LT+SgfPbzQhujyQJDBojETxQ0ERcE713IyMwyc8ysZhXO0TidnhwCARvD59i6ZjJzy6zGzxq/H4Itjbtnwr/v57ef//KrE5a3k5BsAAAs+V43dfnaLmY/6WTau/8Jzbtlf7k5OyvCuTX/jrjPwO34d/9otSU8X7EWkdgzsRzLCXKnxEpFE6nmJSO105m1U4yUitVPubNhpqPESkUU6U59ptFFEasaMOm0sVdTrjSbZiMIvnPeTs8EqZ/0PQZTBYm5NkHIiqr4zmcb0tf4wPoOQjaET/i/w/PXTbpm9lRxiMTvn72vuKr+Og2f82Ia3/3qNW4brkg/yqqP++7z+J34dJz/iv2mrXvG/UrzlPW7ZzV9Onkzj//3TLe468//WP1ZYNH/P8tU+SJXkFgB/AuAaAG0Au8zsP3bLjgI4g07003yKGUREpAY6+bzqf81rHsDvm9mPSa4E8DzJJ8zs5W757WZ2on9VFJHi9SeTap6WbLzM7DiA493/nyF5EMAmAC+HK4pIbXVCJerf8/o5klsB3ALg2e4iA/A4SQPwx2a2y1lvJ7pn4c21a7LWVUQKclnd20hyHMB3AXzWzE53F99mZpMkNwJ4guQhM3v60nW7jdouABi+bkvGmxFFpEhVT4mTqnYkB9FpuL5hZt+7uNzMJrv/TgHYA+DWflRSRIrVSYnDVI+ypBltJICHABw0s68sWD4GoNG9DjYG4A4AX+ypNtFxyNJfa/sbDCfZCLINeGEIrTF/gwzq0TwX/H4M+9tsvuW/dc0LyfsLX3NwhtBwXjMADLw24pbNrUreYXPa315jzB/+Hz3hv+bz7wrWezV5EosL6/0P1du/4pc1z/vv2Yo3/IN89HP+erNTGxKX/9+7/oO7zt2/+y/dsjyySlwO17xuA/DbAP6C5AvdZf8GwCEAezptGwYAfNPMHu1HJUWkWJ2sEtU+bUwz2vgM/D7R+/KtjohUQef2oJo3XiJyJboMel4icmW6HCLsReQKc3G0scrUeIlIIp02LkcUDuH9CITr+IXhTObR5BzO5BbRRBrWCCoZhC803vYrGf0ozq9w9hesM3wyGMZf41dy4Ly/URtLzophTT8uY3zcz1JxYf2oW8Zz/rGaG08+HquP+O9Le1WQ3WLFnFsG+HUcfXqlWzY2m1yXD669z12nebs/MQqe8IvSUA57EaklAzCvnpeI1JFOG0WkfkynjSJSQ3VIRljtfqGIlKbd7X0t9VgKyS0knyR5kOQBkvfnUb/69LyKTKQTpI73mvtoRLE57f9GtEb99RjUwwaCkVTnhu7mGX+Ub/pqf2feCCsADASjfCOvJ4+GRVMEnD2y2i0bSr6/urPN4Obx1jXJEwycP+5vcPSIP5I3t9K/c3/qV4P3Zcg/xsNXn09e3vBHes9v6N/XN+dkhEtlY86kPo2XiBTGQMy38zkx61c2ZjVeIpJoGde81pPcv+D5riCr8la8MxtzZmq8RGQxW9Zp44k0M4c52ZgzU+MlIovkPQGHl425F2q8RCRRXo2Xl425V2q8RGQRA9HK6YI9nGzMZra3l41Wq/HK0tAH+eHjfQXj9dEvTjt5Pc4HNyhHVWwG9WgF2xwKQjNOJ7+tjWB2+HYQstEIQj3aw8HNzV44R3RAgu/L/EgQVnK1f0P30KEVicsHpv3tzazz6zi/xg95GHnD/0rNbJ3xt3lk3C3zcEVwV38O8gpSXSIbc2bVarxEpBJseRfsS6HGS0QSmRovEakf3ZgtIjWlnpeI1I4Z0Mo6GFYQNV4ikqjqKXHq33hFIQ/hekFZlHPeKwqG+NsD/pB2FGIR1XHwLT9DxPx48v5aQ349oowTrbFgvQv+el5IRNvJegEADH7t20F4CKZG3KKZDcmhDa2RKNuHX8fRCf9rMz7h13Fkys9iMXg+eb23d5xz1+Hx5BCQPBiqf9q4ZBRalIuH5A6Sr5A8TPKB/lZVRIqTLpdXmRf10/S8EnPxAHgFwIMAPg5gAsA+ko/0mqNHRKrBisyhl8GSjVeQi2c1gMNmdgQASD4M4C70mKNHRKqh9qeNC12Si2cTgNcXFE90lyWtt5PkfpL7W2fPZqyqiBSlM9rYSPUoS+o9J+TiSWqWEzuaZrbLzLab2fbm+PLv4RKR4pmle5Ql1Wijk4tnAsCWBX+2GcBkvtUTkbJU/bRxycYryMWzD8BNJG8AcAzA3QDu7ak2Rbbi0Q35YRhF8mILs1QERdkiJdxwCAAYPpkcvjC3MuPOgtCRubV+hgUvDGTolN/h98IaAGD4TT8so/kePzHnzNGVicsHz2Q7Hq0gk8abH/RTdwyt9rNKtI+OJS4f3eefqcys698XxsDKN15pThsv5uL5CMkXuo87zWwewH0AHgNwEMC3zexAH+sqIgWylI+ypBltdHPxdJOJ9ZRQTEQqyADT7UEiUkdVP21U4yUiiWofpCoiV5463NuoxktEFjMsMflC+erTeHnHMWvXNkM4BABwLnlFizJRBIkXOOeXmTeBRVAPAJhZmxxG0ZwOJvSIJgIJjlXztP/i6EQ9ROEQI8f9j2Q0gYi9sNovW5d8PGZXB8c3CqUJPh8jk4N+4TG/rL0iuS4rP/6Gu07zkWv8feVAp40iUkPUaKOI1JR6XiJSO6YL9iJSV+p5iUg9qeeVXpHHKvpV8QfD/FG54OJmOKIYjFJGo42NmWDkcDB5PQve7ShPfeNcMFwa1L81mlw25Nw4DgBzq/x6WDDKN3zSL2w4I7PRjdkXNvlDmwPBCOvsWn+94RPByOx15xOX29c2uuuce3+fu0bRiGsFVKvxEpFqUJyXiNSV4rxEpJ7UeIlILem0UUTqKOt8zkVR4yUiixnDEfQqUOOVJMN7FqawD8IJ2PJ31jwdTEc/Eoxje9EcQ349mmeC/PB+6vUwfKHt5Hqf3eiHE4Q3egfHeOZXLvjrTQ0nLr+wJci/PxLEywR1jMIh5sb9F9A+n3zT9s9u8A/w2M2n3LJcVLznVd6kayJSbTkmsSe5g+QrJA+TfCCP6qnxEpFkOTVeJJsAHgTwCQDbANxDcluv1VPjJSKLXQxSTfNY2q0ADpvZETObBfAwgLt6raIaLxFJREv3ALCe5P4Fj52XbGoTgNcXPJ/oLuuJLtiLSLL0F+xPmNn2oDype9bzcIAaLxFJlGOc1wSALQuebwYw2etG1XglyfCmeZkcgDjffLSv1oogHCIKsbiQXBaFSrSD0ItWMMSf5Vjxgh9OYFECiyCTRvutIX9/zvLRY/7Hf27cr0j0vrAdXInZ4MecNN8YSVx+4eZpdx17bp2/rzzkF2G/D8BNJG8AcAzA3QDu7XWjqRovkrsBfBLAlJm9d8HyowDOoJNEZn6JrqOI1MUywiCW3JTZPMn7ADyGznQ0u83sQK/bTdvz+hqA/wTgTxLKbjezE71WREQqJscgVTPbC2BvfltM2XiZ2dMkt2bdSXf0YScANNeuyboZESlQOP1bBfQaKmEAHif5fMLw6C/+yGyXmW03s+3N8fEedykihcgxwr4fer1gf5uZTZLcCOAJkofM7Ok8KiYi5VkQw1VZPfW8zGyy++8UgD3oRNKKyOUgvwj7vsjc8yI5BqBhZme6/78DwBdzq9mlsvwKRMc171+V6PpAFGkQTLIRbTP6VWyNJq/IIMVJY9Yva0VhIPPBNr3QhuADH4VzRCxYr+Uc47kghKU14m9v4Kz/m+9l0gCA4UOjbtns2uT3zIKQmNaKPneNKt7zShsq8S0AH0bnNoAJAF8A8CSAPSQvbuebZvZon+opIgWr+mlj2tHGe5yi9+VYFxGpCqv+aKMi7EUk2eXQ8xKRK5AaLxGpo6pf81I+LxGpJfW8lssb5s8aDpHx1y0Mr3EK20P+FdiBs34Whca0/xtnzSBThRO+0JhzVwGCyUoseY6Kzmrngzo6oR7R++Jl5gDicIjWuD9xRzsKsRhN3ubokeTJQwBgfkyhEiIi76TRRhGpLfW8RKRuiOpfsFfjJSLJ1HiJSO3UIKtEfRovb/AnOsBRWeab4bMkbc+4uWi94JPl3SzNmWC0K7ixOczPH9zQ7Y3mRXnqo9FLzvn1ZzRNgLPJ6Eb1MIgoeM8aQX7+VjBKOXQieb3pm/wc9gxGZnOhC/YiUkfqeYlIPanxEpHaKTnFcxpqvEQkkU4bRaSe1HiJSB3p9qCq6ksYhaMfH4JgmD+8EdzBIFc6/HuNYRlCCqIbxCNZwzm8sijMA0E4B73c/ED8uQqO1dxqJ4f9tF+RoTf6+PXVNS8RqSMi/9/wvKnxEpFk6nmJSB1ptFFE6kmNl4jUjpIRikhtqeeVE+9AZs3YkGVfWUUJ56MLCxmzSrirBOEQWcIrAADtKD+/k91iLqjHoP9zb0FWhug9GzjthBsE60R56gfO+/WfHw9CNqJ8/6vmE5cPHfcT98+u7m/rUvVrXqlmDyK5m+QUyZcuWb6D5CskD5N8oD9VFJFSWMpHSdJOffY1ADsWLiDZBPAggE8A2AbgHpLbcq2diJSGlu5RllSNl5k9DeDUJYtvBXDYzI6Y2SyAhwHclbQ+yZ0k95Pc3zp7tqcKi0gBDJ07Q9I8ekRyC8knSR4keYDk/WnW62XS2U0AXl/wfKK7bBEz22Vm281se3N8vIddikgRLk7AUVDPax7A75vZzQB+DcA/T3MW10vjlXTVsuKX+EQktYKueZnZcTP7cff/ZwAchNMRWqiX0cYJAFsWPN8MYLKH7YlIhdBL/r/YepL7FzzfZWa7Mu2T3ArgFgDPLvW3vTRe+wDcRPIGAMcA3A3g3h62F8syAUfe+8q6v2iihGh70fWEqI5OxokwYsOZtAOIh/jzFtUjKouOR2s0+UA2oslDgtc8v8LfV3M62GZwntP4WfJXcfYqP6WHGwKSh+X1qk6Y2fboD0j+GYBrEoo+b2Z/2v2bcQDfBfBZMzu91E5TNV4kvwXgw+i0sBMAvmBmD5G8D8Bj6CQQ2W1mB9JsT0SqL8+RRDP7WLgvchCdhusbZva9NNtM1XiZ2T3O8r0A9qbZhojUS1G3B5EkgIcAHDSzr6Rdr5cL9iJyOSsuSPU2AL8N4CMkX+g+7lxqpfrcHiQixSkwANXMnkGG3IdqvEQkWcUDn9R4icgiF4NUq6z+jVfWRNtFvjFZ95U1ZMP71EXbiybgCCb7QJiNYvkTX1g0+h9coQ3DQJzXljl0JJgIpJXxKrJXl+ZZf4PRhCR5YJQxpALq33iJSP40e5CI1JUyqYpIPannJSJ1pAv2IlI/BiD9jdmlUOMlIol0zatMYThBzvsqeiKQLPuLkltEmSOievhJDzLti1FYRlSPIGTDnFAJ9mE++yj8IjwPcyIi2iV9QxXnJSL1ZKbTRhGpJ/W8RKSe1HiJSB2p5yUi9WMAWtVuvS7vxqsfI4BFba+X/XlDXtFPaUVGZsNrxNFN4BlG+cKbwDPexB6FF3g3iAOAOd9EC15zI8iXnwf1vESknjTaKCJ1pJ6XiNSPUuKISB0RAHXBXkTqaBkzZpdCjZeILKbTxgrLe5S5H290kTnsl5+Kfmne/qJsBdFV4ijPfiMKA8nyZkdhGdniSrLknLfgdbWHlr255ez58h9tJHkUwBl08gvMm9n2XrcpIuW7UkYbbzezEzltS0Sq4HLveaVBcieAnQDQXLumiF2KSC+s+qONGWeZewcD8DjJ57uN1OI/MNtlZtvNbHtzfDyHXYpI31nKR0ny6HndZmaTJDcCeILkITN7OoftikiJqh4q0XPPy8wmu/9OAdgD4NZetykiFXAxm+pSj5L01PMiOQagYWZnuv+/A8AXc6mZFJrDvtAMHEXzwg2ifPkZj0eUBSLL8WcU5pEpBCQlQxzSUgG9njZeDWAPOzMZDAD4ppk92nOtRKRUhFX+tLGnxsvMjgB4X051EZEqaVe763XlRtiLiO8KOG0UkctU1U8b84jzEpHLUcGjjSSbJH9C8vtp/l49LxFJUEoYxP0ADgJYleaP1XiVrSohCkWGUWTdXpg5Iks9Mh7gfrxn3npZM3D0quDZg0huBvCbAL4E4F+lWUeNl4gkWsY1r/Uk9y94vsvMdi1zd18F8AcAVqZdQY2XiCRL33idWCoVFsk/A3BNQtHn0UmnNWVmz5P8cNqdqvESkcUMQDu/00Yz+5hXRvLLAD5F8k4AIwBWkfy6mf1WtE2NNopIgpQjjTlc1Dezz5nZZjPbCuBuAD9cquEC1PMSEU/F47zUeInIYgagVXyIvZk9BeCpNH+rxisv/ZjcIu+kAX1MQrAseYcTAPWYUCVSlffm5wywat8fpMZLRJLptFFEaifn0cZ+UOMlIsnU8xKRWlLjJSK1Ywa0WmXXIqTGqwhFjyRVbuTqElWpX9E5/avyutNSz0tEakmNl4jUj2m0UURqyABTkKqI1FIJtwcthxovEVnMTFOfiUhN6YK9XNbqNvy/UNa6V+U197keVvGeV8/JCEnuIPkKycMkH8ijUiJStuKSEWbVU+NFsgngQQCfALANwD0kt+VRMREp0cUbs9M8StLraeOtAA6b2REAIPkwgLsAvLzwj0juBLATAJpr1/S4SxHpNwNgFb89qNfTxk0AXl/wfKK77B3MbJeZbTez7c3x8R53KSJ9Z91khGkeJem155V0ybDaQxQikopd5hH2EwC2LHi+GcBkj9sUkSqoeIQ9rYfRApIDAP4SwEcBHAOwD8C9ZnYgWOdNAK8tWLQewInMlciP6vFOVahHFeoA1LMe15vZhqw7Ivlod39pnDCzHVn3lVVPjRcAdCeK/CqAJoDdZvalZa6/f6nZdougelSvHlWog+pRXT0HqZrZXgB7c6iLiEhqmjFbRGqpCo3XrrIr0KV6vFMV6lGFOgCqRyX1fM1LRKQMVeh5iYgsmxovEaklNV4iUkulNV5VSaVD8ijJvyD5Asn9Be97N8kpki9dsrzQYxPUo5BjQ3ILySdJHiR5gOT9C8qKPhZRXYo6HiMknyP5024d/mhBWSW+N5VgZoU/0AlofRXAjQCGAPwUwLaS6nIUwPqS9v0bAD4A4KUyj01SPYo8NgCuBfCB7v9XonPXxraSjkViXQo+HgQw3v3/IIBnAfxalb43VXiU1fP6eSodM5sFcDGVzhXFzJ4GcOqSxYUfG6cehTGz42b24+7/zwA4iE52kjKOhVeXwljH2e7Twe7DoO/NO5TVeKVKpVMQA/A4yee7ecfKdkUfG5JbAdyCTm+j1GNxSV2AAo8HySbJFwBMAXjCzEo/HlVTVg77KqXSuc3MJkluBPAEyUPdnkhZrthjQ3IcwHcBfNbMTpMs7VhcWpfu4sKOh5m1ALyf5BoAe0i+F9X6bJSurJ5XZVLpmNlk998pAHvQ6ZqX6Yo8NiQH0WksvmFm3+suLuVYOHUp5bNiZm8DeArADlTos1EFZTVe+wDcRPIGkkMA7gbwSNGVIDlGcuXF/wO4A8BL8Vp9d8Udm24P6yEAB83sKwuKCj8WXl0KPh4buj0ukBwF8DEAh1CRz0ZllDVSAOBOdEZyXgXw+ZLqcCM6IzY/BXCg6HoA+BaA4wDm0PlV/UwZxyapHkUeGwAfQuf050UAL3Qfd5Z0LBLrUvDx+NsAftKtw0sA/t2CstK/N1V56N5GEaklRdiLSC2p8RKRWlLjJSK1pMZLRGpJjZeI1JIaLxGpJTVeIlJLarxEpJb+P8oc8i1d7QZPAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "im = ax.imshow(np.transpose(r.v()), origin=\"lower\")\n", + "fig.colorbar(im, ax=ax)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "76929fea-527a-4e7e-9bf2-d0af4c320d07", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0.000e+00 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00 0.000e+00]\n", + " [ 0.000e+00 0.000e+00 1.735e-18 ... -4.441e-16 0.000e+00 0.000e+00]\n", + " [ 0.000e+00 1.735e-18 -3.469e-18 ... -8.882e-16 -1.332e-15 0.000e+00]\n", + " ...\n", + " [ 0.000e+00 0.000e+00 6.661e-16 ... 7.461e-14 2.789e-13 0.000e+00]\n", + " [ 0.000e+00 0.000e+00 -1.332e-15 ... 3.464e-14 -2.887e-13 0.000e+00]\n", + " [ 0.000e+00 0.000e+00 0.000e+00 ... 0.000e+00 0.000e+00 0.000e+00]]\n" + ] + } + ], + "source": [ + "print(r)" + ] + }, + { + "cell_type": "markdown", + "id": "5a479f57-aced-49d1-8f27-ae7e34cae5a9", + "metadata": {}, + "source": [ + "## Convergence " + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "e7043f56-2038-4d79-9082-1d7641ca98ee", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "N = 32, error = 8.044915834004946e-05\n", + "N = 64, error = 2.0120304937333502e-05\n", + "N = 128, error = 5.03080680318748e-06\n", + "N = 256, error = 1.2577618783599838e-06\n" + ] + } + ], + "source": [ + "for N in [32, 64, 128, 256]:\n", + " mg = MG.AxisymmetricMG2d(N, N,\n", + " xl_BC_type=\"dirichlet\", xr_BC_type=\"dirichlet\",\n", + " yl_BC_type=\"dirichlet\", yr_BC_type=\"dirichlet\",\n", + " xr_BC=rr_func, yr_BC=zr_func, true_function=true,\n", + " verbose=0)\n", + " mg.init_zeros()\n", + " mg.init_RHS(f(mg.x2d, mg.y2d))\n", + " mg.solve()\n", + " v = mg.get_solution()\n", + " e = v - true(mg.x2d, mg.y2d)\n", + " print(f\"N = {N}, error = {e.norm()}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b21f2c2-c105-44b9-9e1e-b6e3e384f1c0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}