{ "cells": [ { "cell_type": "markdown", "id": "eb519d17", "metadata": {}, "source": [ "\n", "# Alternative Linear Solvers (Scipy and Petsc)\n", "Here we look at different ways of solving PDEs using external\n", "packages and python functionality.\n", "Different linear algebra backends can be accessed by changing setting the\n", "`storage` parameter during construction of the discrete space. All\n", "discrete functions and operators/schemes based on this space will then\n", "use this backend. Available backends are `numpy,istl,petsc`. The default is\n", "`numpy` which uses simple data structures and linear solvers implemented in\n", "the `dune-fem` package. The simplicity of the data structure makes it\n", "possible to use the buffer protocol to seamlessly move between C++ and\n", "Numpy/Scipy data structures on the python side. A degrees of freedom\n", "vector (dof vector) can be retrieved from a discrete function over the\n", "`numpy` space by using the `as_numpy` method. Similar methods are available\n", "for the other storages, i.e., `as_istl,as_petsc`. The same methods are\n", "also available to retrieve the underlying matrix structures of linear\n", "operators.\n", "\n", "We will mostly revisit the nonlinear time dependent problem studied at the end of the\n", "[concepts section](concepts_nb.ipynb) which after discretizing in time had the variational formulation\n", "\\begin{equation}\n", "\\begin{split}\n", "\\int_{\\Omega} \\frac{u^{n+1}-u^n}{\\Delta t} \\varphi\n", "+ \\frac{1}{2}K(\\nabla u^{n+1}) \\nabla u^{n+1} \\cdot \\nabla \\varphi \\\n", "+ \\frac{1}{2}K(\\nabla u^n) \\nabla u^n \\cdot \\nabla \\varphi v\\ dx \\\\\n", "- \\int_{\\Omega} \\frac{1}{2}(f(x,t^n)+f(x,t^n+\\Delta t) \\varphi\\ dx\n", "- \\int_{\\partial \\Omega} \\frac{1}{2}(g(x,t^n)+g(x,t^n+\\Delta t)) v\\ ds\n", "= 0.\n", "\\end{split}\n", "\\end{equation}\n", "on a domain $\\Omega=[0,1]^2$. We choose $f,g$ so that the exact solution\n", "is given by\n", "\\begin{align*}\n", "u(x,t) = e^{-2t}\\left(\\frac{1}{2}(x^2 + y^2) - \\frac{1}{3}(x^3 - y^3)\\right) + 1\n", "\\end{align*}\n", "The following code was described in the [concepts section](concepts_nb.ipynb)" ] }, { "cell_type": "code", "execution_count": 1, "id": "fffb880a", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:06:21.168390Z", "iopub.status.busy": "2024-03-01T11:06:21.167996Z", "iopub.status.idle": "2024-03-01T11:08:02.456303Z", "shell.execute_reply": "2024-03-01T11:08:02.454826Z" } }, "outputs": [], "source": [ "import numpy, sys, io\n", "import matplotlib.pyplot as plt\n", "\n", "from dune.grid import structuredGrid as leafGridView\n", "from dune.fem.space import lagrange as solutionSpace\n", "from dune.fem.scheme import galerkin as solutionScheme\n", "from dune.fem.function import gridFunction\n", "from dune.fem import integrate\n", "from dune.ufl import Constant\n", "from ufl import TestFunction, TrialFunction, SpatialCoordinate, FacetNormal, \\\n", " dx, ds, div, grad, dot, inner, sqrt, exp, sin,\\\n", " conditional\n", "\n", "gridView = leafGridView([0, 0], [1, 1], [4, 4])\n", "space = solutionSpace(gridView, order=2)\n", "\n", "x = SpatialCoordinate(space)\n", "initial = 1/2*(x[0]**2+x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1\n", "exact = lambda t: exp(-2*t)*(initial - 1) + 1\n", "\n", "u_h = space.interpolate(initial, name='u_h')\n", "u_h_n = u_h.copy(name=\"previous\")\n", "\n", "u = TrialFunction(space)\n", "v = TestFunction(space)\n", "dt = Constant(0, name=\"dt\") # time step\n", "t = Constant(0, name=\"t\") # current time\n", "\n", "abs_du = lambda u: sqrt(inner(grad(u), grad(u)))\n", "K = lambda u: 2/(1 + sqrt(1 + 4*abs_du(u)))\n", "a = ( dot((u - u_h_n)/dt, v) \\\n", " + 0.5*dot(K(u)*grad(u), grad(v)) \\\n", " + 0.5*dot(K(u_h_n)*grad(u_h_n), grad(v)) ) * dx\n", "\n", "f = lambda s: -2*exp(-2*s)*(initial - 1) - div( K(exact(s))*grad(exact(s)) )\n", "g = lambda s: K(exact(s))*grad(exact(s))\n", "n = FacetNormal(space)\n", "b = 0.5*(f(t)+f(t+dt))*v*dx + 0.5*dot(g(t)+g(t+dt),n)*v*ds" ] }, { "cell_type": "markdown", "id": "b479a9c1", "metadata": {}, "source": [ "When creating a scheme, it is possible to set the linear solver as well as\n", "parameters for the internal Newton solver and the linear solver\n", "and preconditioning. See a list of available solvers and preconditioning\n", "methods at the end of this section." ] }, { "cell_type": "code", "execution_count": 2, "id": "2ab2c1eb", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:08:02.463821Z", "iopub.status.busy": "2024-03-01T11:08:02.463123Z", "iopub.status.idle": "2024-03-01T11:10:06.166593Z", "shell.execute_reply": "2024-03-01T11:10:06.165111Z" }, "lines_to_next_cell": 1 }, "outputs": [], "source": [ "scheme = solutionScheme(a == b, solver='cg')\n", "\n", "endTime = 0.25\n", "exact_end = exact(endTime)\n", "l2error = gridFunction(name=\"l2error\", expr=dot(u_h - exact_end, u_h - exact_end))\n", "h1error = gridFunction(name=\"h1error\", expr=dot(grad(u_h - exact_end), grad(u_h - exact_end)))" ] }, { "cell_type": "markdown", "id": "64044305", "metadata": {}, "source": [ "We define a function to evolve the solution from time 0 to the end time.\n", "The first argument is a class with a `solve` method that moves the\n", "solution from one time level to the next - i.e., solves the non-linear\n", "problem for $u^{n+1}$ given $u^n$:" ] }, { "cell_type": "code", "execution_count": 3, "id": "0204e0ab", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:10:06.172482Z", "iopub.status.busy": "2024-03-01T11:10:06.172223Z", "iopub.status.idle": "2024-03-01T11:10:06.177819Z", "shell.execute_reply": "2024-03-01T11:10:06.176708Z" }, "lines_to_next_cell": 1 }, "outputs": [], "source": [ "def evolve(scheme, u_h, u_h_n, endTime):\n", " time = 0\n", " while time < (endTime - 1e-6):\n", " t.value = time\n", " u_h_n.assign(u_h)\n", " scheme.solve(target=u_h)\n", " time += scheme.model.dt" ] }, { "cell_type": "markdown", "id": "a3238df5", "metadata": {}, "source": [ "We can simply use the `galerkinScheme` instance `scheme` in this function to produce the solution\n", "at the final time. We combine this with a loop to compute the error over\n", "two grids and estimate the convergence rate:" ] }, { "cell_type": "code", "execution_count": 4, "id": "e9e42c41", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:10:06.182501Z", "iopub.status.busy": "2024-03-01T11:10:06.182258Z", "iopub.status.idle": "2024-03-01T11:10:47.885707Z", "shell.execute_reply": "2024-03-01T11:10:47.884426Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer: step: 0 , size: 16\n", "\t | u_h - u | = 1.55256e-04 , eoc = -\n", "\t | grad(uh - u) | = 3.99906e-03 , eoc = -\n" ] }, { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrASEDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3yV2jhd1ieVlGRGhGW9hkgfmRWTL4ltItOs73ybh47m2N3hVXMcIClnbLdt65AyeeAa1pY1mieNiwVwVJRip/AjkfUVmJ4a0pLOC0EMxggXZGr3MrYTABTJbJQhVyp+U45FAFm/1nS9KKDUdSs7Pf9z7ROse76biM1S/4TDwx/wBDHpH/AIHRf/FVhfE2PzdI0pf+oiP/AERNXnn2L2ryMdmn1Wr7Plvpfc9vL8phi6XtJTtrbb/gnsH/AAmHhj/oY9I/8Dov/iqP+Ex8Mf8AQx6R/wCB0X/xVeO/Yvam/YvauT+3v7n4/wDAO/8A1dpf8/fw/wCCeyf8Jj4Y/wChj0j/AMDov/iqP+Ex8Mf9DHpH/gdF/wDFV4x9i9qb9i9qf9u/3Px/4BX+rlL/AJ+/h/wT2n/hMfC//QyaP/4HRf8AxVJ/wmXhf/oZNH/8Dov/AIqvFfsXtUEdl9/j+M1Szxfy/iH+rdK/8X8P+Ce4/wDCZeF/+hk0f/wOi/8AiqP+Ey8Lf9DLo/8A4HRf/FV4a1l7VE1l7VSzpP7P4mi4YpP/AJffh/wT3f8A4TPwt/0Muj/+B0X/AMVSf8Jn4V/6GXRv/A+L/wCKrwVrL2qJrL/ZqlnC/l/EtcK0n/y+/D/gnv8A/wAJp4V/6GbRv/A+L/4qk/4TXwp/0M2jf+B8X/xVfPbWI/u1WnsAIjxWizVP7I3wnTtdVvw/4J9G/wDCa+FP+hm0b/wPi/8AiqP+E18Kf9DPov8A4Hxf/FV84Np49KibTzzVLM0+g/8AVGPSr+H/AAT6U/4TXwp/0M+i/wDgfF/8VR/wm3hT/oZ9F/8AA+L/AOKr5maxYdqiazYH7tWsxT6EvhL/AKefh/wT6e/4Tbwn/wBDPov/AIHxf/FUf8Jt4T/6GfRf/A+L/wCKr5dNqR/DUZt/aqWPXYzfCjX/AC8/D/gn1L/wm/hP/oaNF/8ABhF/8VR/wm/hP/oaNF/8GEX/AMVXyt9n/e9O1H2f2p/XV2M/9WJfz/h/wT6p/wCE38J/9DRov/gwi/8AiqP+E38J/wDQ0aL/AODCL/4qvlX7P7U37P7U/rq7C/1Zl/P+H/BPqz/hN/Cf/Q0aL/4MIv8A4qj/AITfwl/0NGif+DCL/wCKr5S+z+1N+z+1P64uxH+rcv5/w/4J9X/8Jv4S/wCho0T/AMGEX/xVbFpeWuoWsd1ZXMNzbyDKSwuHRh04I4NfG/2f2r6i+GYx8OdFHpC3/obVtRre0bR5uY5Y8HGMm73N3Utb0nRvK/tTVLKx83Pl/arhIt+MZxuIzjI/MVRXxr4Udgq+JtGZicAC/iJJ/wC+q8w/aAXdP4b/AN26/nDXkOnx41K1P/TZP5ilOvyz5bFYbK/b4f23Nbfp2Pry61WC3tre5RWubeeaOESwMrKu9gisckZGSBxk+1UrzxRY2V/c2kkdwzWyFnkVBs3BN/lgkjLlSDj0+hrQ1DTbbU4EhuhKUSRZQI5njO5TlTlSCcEA46cVXk8PaZMzNLA8hZNjb5nbd8mzccnl9vy7/vY710HkFWPxVZSX0Fmbe6SeRikgZVxAd7IN7BsDcykDBOfxFXNF1iPW7E3KWt1asr7Hgu0CSxnAbDLk44YH8aYvh7TEaNhA+6M7gTM5LHcXBb5vnIYlgWzgnIxVuy0+309JFt1f94/mO0krSMzYAyWYknhQOvagCzRRRQAUUUUAcl49j82z0lf+n/8A9oTVyn2L2rtvFsfmjR1/6fj/AOiJqzDZcdK+D4lrcmNS/ur82e7luI9nR5fM5r7F7U37F7V0psuOlNNlx0rwFij0ljDmvsXHSm/YvaukNlx0ppsuOlUsUUsYc39i9qrQ2WfM4/5aGurNlx0qrBZcTcf8tGrSOK0Zaxmpz7WXtUTWXtXTtZe1RtZe1WsUbRxnmcu1l7VE1l7V1DWXtUTWXtWqxRtHGeZzDWXtVW5sv3Lcf5zXWNZe1Vbuyxbvx6fzraGK1Rr9c0Oday9qiay9q6hrL2qJrL2q1ijeOMOXay9qjay9q6drL2qJrL2rVYo2jjPM5hrL2qFrEc/LXUNZe1RtZe1axxRtHGHJtYDzz8v8NB08Y6V0Zsv9JIx/B/WnGy46Vr9a8yo4pHMHTxjoaYdP4rpzZcdKabLjpVLFvuafWYnMHTzimnTziunNlx0ppsuOlUsWP28H0OXNg2K+i/hwu34f6QvpG4/8favGjZcdK9p8ALt8EacvoJB/5EavXyqt7SckfK8VTjKjT5e7/I89+PK7rjw5/uXX84a8oso8X1ucf8tV/nXr/wAb4/Mu/Dw9I7r+cNeXW9vsuYmJwFcEk/WtMTO1e3odOTUVLKub/F+p9YUVRutVgt7a3uUVrm3nmjhEsDKyrvYIrHJGRkgcZPtUEuuxw6rNp7WV4HSIypKUVYpsAEqjMwBbn2xg16x8CatFc/B4w0+eWzj8qdPtfEbMY8ZLsgAIY7ssvVdwwVJIBBrV02/XUrMXKwyw/vJI2jl27lZHZGB2kjqp6GgC3RRRQAUUyWWOCJ5ZpFjjRSzu5wFA6kk9BVd9V06OGCZ7+1WKc7YXMyhZD6Kc8/hQBmeJTiTRv+v4/wDoiaoNwqTxUcHRz/0/H/0RNVPzK+J4hwntsWpeS/NmFTFeylylncPQUZU9qreZS+ZXgPLn2Esw8yf5fSjap7VB5lL5lQ8vZoswfcm2J6VWtokPnf8AXVqk8yoLWT/Xf9dWqfqMkmaRzB9ywbdT6U02intTxJS76xeFqI3jmLK7WY9Kiay9qvb6XeKh0qqOiOYmY1l7VTvbL/RZOPT+YroMg9qr3oU2knHp/MU4uopI3jmCMlrL2qJrLnpXQmNDTTAh9KXtpLc6I49dzm2suOlRtZe1dKbVTUbWY9KpYlnRHHLucy1l7VE1l7V07WXtUbWXtWixRvHG+ZyTWX+mdP4P6082XHSt42X+nYx/yy/rUpsvatnitjSOMOaNl7U02XHSukNl7U02XtTWKNFjDmzZe1MNlx0rpTZe1NNl7VSxRSxnmc0bL2r0vwOu3whZL6NKP/Ir1y7WXtXWeD12+GbdfSWYf+RXr6bhut7SpNeSPIzmv7SnBeZxHxgTff6CP+mV1/OGvOkgIdT6EV6b8VE8zUtDGOkN1/OGuEaDCk46CurH1LYy3oe5k9Xlyvl/xfqe+ahpttqcCQ3QlKJIsoEczxncpypypBOCAcdOKiOjWbSrK32lnWPy1ZrqU7RjGR83DYH3hzyeeTWhRX0h8GZ8WiWEPkbI5cQO0iK08jLvLFizAthmyxOTk5qXT9NttLt2gtRKI2kaU+ZM8h3McscuSeSSfxNW6KACiiigCpqVl/aGny2vmeWXwVfbkBgQRkdxkDI71g3XhK4u7SOF9SRSks8hZIGUHzW3sMCTnDbuGypBAZWxmupooA5jxejCbRpfOcKLwr5WF2k+RLz0zn8cVn+ZV/xsXEGkGPbu+3/xdP8AUTVgb7r/AKY/rXFiMJ7aXMfLZzifZ4lRv0X6mh5lL5lZ2+6/6Y/rS77r/pj+tcjy1djy1jn3NHzKPMrP33X/AEx/Wjfdf9Mf1rN5ajRY59zR8yobWT/Xf9dWqrvuv+mP61DbPdfvf9T/AKw+tQ8tVjSOOfc2BJTvMrND3f8A0x/WlD3f/TH9a5pZabRx3maXmUokrND3f/TH9ad5l3/0x/WueWWnRHHeZpCSobyT/RH/AA/mKqB7v/pj+tRXb3X2V/8AU9vX1rneW6nRDHeZseZThJWaHu/+mP60u+7/AOmP61zSy03jjvM0g9KHrN33f/TH9advu/8Aph+tc0stN447zNHfS7h6Cs7fd/8ATH9aUPd/9MP1rnllxvHHvuWML/aHT/ll/Wp9qHtWVvu/t/8Ayx/1Xv61Y33frB+tZyy9m6x77lzYlHlIaqb7v/ph+tG+7/6YfrWTwEu5oswfctGFKQ26n0qvvu/+mH60u+7/AOmH61LwUu5azB9yY2y1oeFht0GMek9wP/Iz1k77v/ph+ta3hXJ8Px7sbvPnzjp/rnr6jhajKnVqXfRfmE8S6yt2OV+JCb9V0Uf9MLn/ANChrjbiEJbSseAEJyfpXeePE36xo4/6d7n/ANChrlryDFlOfSNv5VeaVLZjb/D+h9Nl1XlwPL6nqN9qgt7O2u7VYbqCa4hiLrNgBZJFQMuAQ2Cw44+tVJvEcY8SnQ7dIJLkW7ykyT7PnGwhMYJ5V9xbt7841Lyws9QjWO9tILlEYOqzRhwGHQgHv70os7UIiC2hCRxmJF8sYVDjKj0XgcdOBX2R8qZOn67dXo0xpLGGNL+KSRGW4LYCnKH7gyGUg54IzjBq5o+oXGowXElxbRQ+XcPCpimMivs4Y5KqR8wZcY/hz3qw+nWMk8E72Vu01uMQyGJS0Q9FOMj8KmiijgjEcMaRxr0VFAA/AUAPooooAKKKKAOW8bHEGkH/AKf/AP2hNWDvrY+ID7LHSW/6iH/tCauUFyfWvSwmGdWnzeZ+f8Tyaxqt/KvzZp7qXdWaLo+tPF0fWt3gX2PnfayL+6l3VRF1Thcisngn2KVdl3dUNs3+t/66GohcLTLaZT5vP/LQ1jLBvsaRxDszQ3Uoaq4cHvTg1c8sIaRxJYDUoaoN1KGrnlhDeOJJ91RXbf6K/wCH86QNUV23+iv+H8655YQ6aeJ1ReDUoaoA1KGrmlhDaOJLAelDVBupd1c08IdEcUWA1LuquGpwauaWDOiOKDd/p/8A2y/rVjdVHf8A6f8A9sv61Y31hLCG/wBaJt1LuqDfS7qyeE8i1iifdRuqHdRurN4QtYon3Vq+Ff8AkAR/9d5//Rz1ibq2vCn/ACL0X/Xaf/0c9etk9H2c5PyPSwFXnk0ZPjBN+t6UMdLa5/8AQoa5/UIMaZdHHSF/5Gup8SJ5mvaYP+nW5/8AQoayNWh26NfNzxbyHp/smvns5qWzW3+H9D6zC1eXD8vqd9RWdfam1tZ211Db+ZFLcQwuJS0TKJJFQHaVzkFhwcfWs3UvF1tYa1NpixxyyxWrzNmYKd42YTbjOCHBLdAPXnH6AeMdHRXNJ4viOp6fpzwwpd3M0sMim4GIynmgFcjL7miI6DHfnAOno+pyanFcma3EEtvOYWQMx5Cq2fmVWHDDqB6jIIJANKiiigAoqlq8l1FpVw9mG88L8pRdzAZ5KrzkgZIGDkjoa5u41LX/ALDam3F3uSWcSM1kfMm2uPJQjbhA6Hl8AKRzg5FAEXxLfZpWlH/qIj/0RNXBi5967L4rySppOk4RPJ+3jL7zuDeTLxjGMe+fwrzYXPvX1OTU+bDN+b/Q+K4hpc+LT8l+bNoXPvThc+9YwufenC5969R0D594c2Rc+9PFz71ii596eLn3qXQJeHNkXPvTba4/1nP/AC0NZQufekt7n/Wc/wAZrOVAn6vozoFufepluj61hLc+9TLc+9YywsX0MXQaN1br1qRbhTWGtz71Ktz71zTwUWZuEkbYlB71HdPm2fn0/nWatz70XFz/AKM3Pp/OuWeB7DhKSkja3Uoas9br3qVbkHrXLPBPsUq0kXd1LuqqJlPenhwe9cs8J5GscSWA1ODVXDU7dXPLCG8cSG7/AE7/ALZ/1qxuqlu/03/tn/Wp91YPCG/1km3Uu6oN1LurN4QpYkm3Uu6oN1LurN4Q0WJJ91b/AIS/5F2H/rtP/wCjnrmt1dL4R/5FuD/rrP8A+jnrShR9m2z6HIavPUmvIh1sA+INNz/z63P/AKHDWfrIX+w9Q/69pP8A0E1f1048Qab/ANetz/6HDWbrDf8AEkv/APr2k/8AQTXxOc0ebNeb/D+h9Qq/KuU6y8sLPUI1jvbSC5RGDqs0YcBh0IB7+9BsbQxLEbWAxpEYVTyxhYzgFAP7pwOOnAqxRX35JXWxs0hihS1gWKIkxoIwFQkEHA7cMw/E+tOtbO2sovKtLaG3jzu2RIEGfXAqaigAooooAKKKKAPO/jC+zw3pjf8AUSX/ANEzV5GLn3r1b42vs8J6Y3/UTX/0TNXiQufevuOHafNg2/7z/JHzua0ueun5GyLn3p4ufesUXPvTxc+9e46B5Tw5si596cLn3rGFz704XPvUugQ8ObQufekgufv8/wAZrJFz70Q3H3uf4jWcqGpLw+jN9bn3qZbn3rCW596lW596h0DCWHN1bn3qVbn3rCW596mW596xdAwlhzcW596Jrj9w3Pp/Oshbn3pZbnMLc+lYzo6My+r6o6Bbn3qVbn3rDW596mW596zlQOeWHNxbn3qVbn3rDW596lW596xlh0+hjLDs3Vuj61Kt161hLc+9Src+9c08HFmTpSRrCdTe/wDbP+tWBMvrWALn/S+v/LP+tWBc+9c7wSCSkrGz5o9aXePWsgXPvThc+9ZvAk3kjW3+9LurLFz704XR9azeBY+eSNPdXVeEP+Rag/66z/8Ao164UXR9a7jwad3ha2PrJMf/ACK9cGLw7pJM+s4Um5Val+y/Mg8QHGv6Z/163P8A6FDWXq740S/P/TvJ0/3TWl4kONe0z/r1uf8A0KGsjVW/4k97/wBe8n/oJr4rMKHNjub0PpMRW5a/L6HS6zrS6dp9perPBBHLdQxMLtChZXkVSACVKsAS3IPTpVa78QNHrMVtBe2LWsts0okIDeXhdysSJOVPXlVXBHz5IB6KivqD0zi7PxbqNxPYI0VntmYq55DTEzPGPKAZgdoVWcgsAGyM451/C+sXesWdw92kO+GURiSAgo/yKxwQ7A4JIznt0ByK3aKACiiigAooooA8u+O/HgzTc/8AQUT/ANEzV4GCvp+te8/Hw7fBOnH/AKiif+iZq+fBJX6HwtTUsC219p/kjzMZG9T5FwFf8mngr6frVMSU4SV9E6MexxuDLgK+n604FfT9apiT3pwkqXRj2IcGXRs9P1oi2/Nx/Ee9VRLSxS/e+tZyoxutCXB2NFdnp+tSLs9P1qistTLJUOkuxhKDLq7PT9amUJ6frVFZalWSspUl2MZQZeUJ6frSyBPKPH61VWWnPJ+6PNYTpKz0MeR3NFRH6frUqrH6fqaorJUqyVm6S7GEoMvKsf8Ad/U1KqR/3f1NUllqVZfesnSXYwlBl1Ui/u/qalVIv7v6mqSy1MsnvWTpLsYyiycJF9p+7/B6n1qcJF/d/U1REn+k/wDAP61OJPesvZLsZyiyyI4v7v6mnCOL+7+pqsJKeJKTpLsZOMu5ZEcX939TThHF/d/U1WEnvThJUOkuxDjLuWRHF/d/U16f4HAHhCyA6bpf/Rr15UJK9U8Dc+DrE+8v/oxq8XOIKMI+p9LwwmqtS/ZFfxOca7pf/Xtc/wDoUNY2qN/xKb3/AK4P/wCgmtbxWca5pf8A17XP/oUNYept/wASm84J/cPx/wABNfI1qHNV5jpzCty47l9D0misbWdYfTNPtLt3gtPMuoYnju8ZKvIqkAh8BgCWzyOKztT8UvbXl0mn3NheLFZG58lThkyAUJbcdwIy2AvC855Ge8+pOqorjovFGotLp4YWEi3DGNvKI+YrI6M65fPlgKGztbIznbWx4Z1O71TTpHvXs5biGXynlsSTA/yq2UJJJHzYPuCO1AGzRRRQBU1K9/s/T5boR+YyYCpu2gsSAMnsMkZPasa88Xw2FjZy3EUST3F4bUxtOFVds3lO4YgZAJBAxk5HTkjoZYo54niljWSN1KujjIYHqCO4qGLTrKCAww2dvHEXVyiRKF3DGDgDqNq4PsPSgDy79oGZR4N02HD7jqSNnYduPKlH3sYzz0zmvnjca+iv2gf+RF07/sKJ/wCiZa+c81+kcJtfUH/if5I46694eHNOEhqLNLmvqNGYWRMJaeJfeq1GaXKieRFsSU6OTr9ap7jSo5GfrUShqiXTNFZKmWWs5ZalWSocDGVM0Vl96lWSs5ZPepVlrJwMZUzRWSnvJ+7NUVlpzSfuzWFSHusxdPU1FlqVZfes5ZKlWWocDCVM0VkqVZKzlkqVZaxcDCVM0VlqZZKzlkqVZKycDCVMuCT/AEj/AIB/WpxLWaJP3/8AwH+tTCWsuQzlTL4lpwk96oCWniWpcDN0y8JKeJKoCX3p4lqXAh0y8JPevX/AZz4LsD/10/8ARjV4oJa9q8AnPgnTj7Sf+jGr5/Po2pw9T3+Ho8tSfoip4uONb0r/AK9rn/0KGsHUW/4ll3/1xf8A9BNbXjRwmtaST/z73P8A6FDXOX8wOnXI9Yn/AJGvDp4fnhzHl5zVcc15f8P6Hq1FFFcp98FFFFABRRRQAUUUUAeT/tBf8iLp3/YUj/8ARM1fOOa+jf2g/wDkRNP/AOwpH/6Klr5wzX3vDNTlwbX95/kjlrL3h+aM0zNLmvpVWMbD80ZpmaXNaKsKw/NAPX60zNCnr9ar2uqCxJmnByKjzS5rdSTJsTrLUyyVTpQxFDimS4JmgslPMnyGqCyetSGTKnmsakPdZk6eppLLUqy1nLJUqye9Q4HPKmaKye9TLLWcstSrLWLgYypmisnvUqy1nLJUqy1k4GEqZeEn77r/AA/1qUSVnCT97/wGpRJ71ioGcqZfElOEnvVES04S0nAzdMvCSnCWqQl96cJPepcCHTLwlr3T4fHPgbTD/syf+jGr5/Ele/8Aw8OfAeln/Zf/ANGNXzXEcbU4erPYyWPLUn6GT8QZPL1bRz6wXP8A6FDXI3VxmznABP7tuPwrpPidJ5ep6Ic9Ybr/ANChriZLjdE4z1UiubAU74O/qeNnFLmzTm/w/oeyaze3UGn2lzAl7DKbqEPDFAJmKGRQ4YKHwNm45BHQc9qz9R1LUmv7gaa2obFszJtksiE3EAjyyYwS4GSQxPJC4Jzt6mivmz7s5ODUNTa7sYpJtUAJmaV2sPkaENIIycR5EhG04yBgZI5AbW8O3F1daaZrprolpGMYu4PKlVOMBxtUZ78DHOMnGTrUUAFFFFABRUVzcw2dtJcTvsijGWOCfyA5J9h1qm2uWCWsVyZJfJlcoHFvIQrBtpD/AC/JhuDux0PpQB5t+0KceA9P/wCwrH/6Klr5u3exr6R/aF/5ELT/APsKx/8AoqWvm2vschk1hXbu/wBDnq/ELu9jS7vY0lFe8pS7mQu72NLu9jSUVopS7gLu9jQG68Gigd/rVKUrrUQ7d7Gl3exptOrqhKXcQu72NLu9jSUtdEZS7ki7vY0Fjt6Gihvu1rK7i9REiykdjUizexqCnA03F9yXFFpZvZvyqVZvZvyqspqVTWbg+5jKKLKz+zflUyz/AOy35VVU1KprJwfcwlFE4n/efdb7vpUon/2W/Kqyn95/wGpQaxUHrqZSiiYT/wCy35U4XH+y35VCDTwaTg+5DiiUXH+y35U4XH+y/wCVRA04GpcH3IcUTC4/2X/Kvov4cnPw/wBJPrG//obV84A19H/Dj/kn+kf9c2/9DavluJ1alT9X+R6WVpKcvQ5j4tyeXqGhH1iuv5w158LjcQPXiu3+NEnl3vh8+sd1/OGvM4rgmZAOSWGM08qp3y6/+L9TgzClzY7m9D6morG1m7vbXT7SdBPHP9qhWSK0iNwGQyKHB+Qnbs3HOF6daimv5zrLokmoLaS2BkXFk2I3yMYPlk7sZ+U56fd7V8YfUm9RXG2+o+JSdI32l48JmmFzKIo8uv73YGDBCoG2M52ruyOnQ6fhG61W60uVtX87zxKApli8slfLQngomfmLfw+2TjJAN+iiigCvfWcd/ZyW0rMquB8yH5lIOQRnuCAfwrM/4RxPKgi/tO+8uG5a62HysSOzbjuGzkbiWA7E8dFxt0UAeQ/tDRKfBOnTZfcNSRcbztx5Up+7nGeOuM184V9J/tDf8iFp/wD2FY//AEVLXzZX1uRu2Gfq/wBDCruLRSUte4pGQtFJS1opCFoHekoHf61alqgH0tNpa6oSJHUtNpa6IyEOoPSkpT0ro5vdYh1KKSitySRTUymqwNSqaloiSLKmpVNVlNTKayaMZImB+f8ACpA1Vwf3n4VKDWKW5lJEwNODVCDTwaTRm0TA04GoQaeDUNENEoNfSPw3/wCSfaR/1zb/ANDavmsGvpT4bf8AJPdH/wCubf8AobV8pxT/AAqfq/yPQy1e/I4X46yeXc+HTnql1/OGvKrS43XsAz1kUfrXpf7QUnlzeGz6rdf+0a8g06fOp2g9Zk/mK7clp3ym/lL9SsTS5q/N6H2bRRRX56esFFFFABRRRQAUUUUAeSftDf8AIhaf/wBhWP8A9FS182V9J/tD/wDIhaf/ANhWP/0VNXzXX1WSu2Hfq/0MKm47NFJRXspmY6ikorRSEOpB3+tFC9/rWilqgHUtJS10wkSOpabS10wkIdQelJQelb83usQ+lptLXVFki09TTKUVYmTqalU1XVqlU1m0ZSRMD8/4VIDUAPz/AIVIDWKW5lJEwanA1CDTwaTRm0TBqcDUINPDVLRDRMGr6X+G3/JPNH/65N/6G1fMgNfTXw1/5J3o3/XJv/Q2r5Liv+FT9X+R35cvekeb/tFHEnhn6Xf/ALRrxnTGP9q2eOT56Yyf9oV7J+0YcSeGfpd/+0a8Y0pv+JvZf9d4/wD0IV15LW5cp5f8X6nVUX7w+w9ZXUJNPtHgtZ2u0uoXdLS5CgIJFL5LFAy7AwwR36d6zNTtfEFzcXUtkLqCSSNjEzTrsRDBgR7AxHmCX5t2MY/ixxXWUV+enWcva6XqhXTPMa9VreWaR2mvWy6At5aSBWIYkspJwRhCO4FXvDlvq1tBeLq7BpWuS8bCcyAqUXO3IG1d27C9q2qKACiiigAoqlq8d1LpVxHZlhOV+UI21iM8hW4wSMgHIwT1Fc7Naa9Jp0Uccd8joblY83a7lZmBgd2D/OiKSrAliSOjcGgDk/2h/wDkQdP/AOwrH/6Kmr5rr6++IXgqTx1p1jpj3otrSK6E8rLHufIjdRgk4xlgMY/GuE/4Zy03/oYbv/vwv+Ne1l+PpUKXJPe5nOLb0Pnylr6C/wCGc9N/6GG7/wC/C/40f8M56b/0MN3/AN+F/wAa9FZxhvP7iPZs+faWvoL/AIZz03/oYbv/AL8L/jUc/wCzrZJbytBrtzJKEJRGhUBmxwCc8c1SznDef3B7NngNA719Bf8ADOmm/wDQw3f/AIDr/jR/wzppo/5mG7/8B1/xq1neFvu/uF7OR8/0tfQH/DOunf8AQw3f/gOv+NH/AAzrp3/Qw3X/AIDr/jW0c+wi6v7heykeAUte/wD/AAztp3/QwXX/AIDr/jUcf7PFmXlEmu3Kqr4jIhU7l2jk88c5H4VtHiHBrq/uF7KR4LSnpXvv/DO+nf8AQwXX/gOv+NH/AAzvp5/5mC6/8B1/xrX/AFkwVrXf3C9jI8Epa97/AOGedP8A+hguv/Adf8aX/hnnT/8AoYLr/wAB1/xrePFGAXV/cL2EjwSivdz+z3aC4RRrtwYijFn8lchsjAxnuC35e9Sf8M9af/0MF1/4Dr/jWv8ArVl/d/cL2EzwYGpVNe6f8M9WH/QwXX/gOv8AjTh+z5YD/mYLn/wHX/Gk+Kcu7v7hPDzZ4YD8/wCFSA17h/wz9Y5z/wAJBc/+A6/402b4BW6oDDrk7tvUENAo+UsNx69hk/hWa4ny/u/uIeFmeKA04Gvbv+FBWX/QfuP/AAHX/Gj/AIUHZ/8AQfuP/Adf8aP9Z8v7v7iHhKh4mDTga9r/AOFCWf8A0H7j/wAB1/8AiqX/AIULaf8AQfuP/Adf/iqn/WbAd39xLwdQ8WBr6d+Gn/JO9G/65N/6G1cZ/wAKHtP+g/P/AOAy/wDxVei+FdHk8P8Ahu10qSQSG2MiK4GNy72Kkj12kZr5/Ps1w+NpwjRvo+qOnC0J0m3I8h/aPOH8MfS7/wDaNeLaS3/E5sf+viP/ANCFfVXxB+HFt8QDpxuNQms/sPm7fLjDbt+zOc+mwfnXHW/7PNhbXMU8fiC5LxOHXdbqRkHIzzWGDzOnRwvsZb6/idEoNyuez0VjazFqB0+0MEc9zdx3ULubSQQgoJFL5DOAVKbhtJPWszUbXXrm6uZbVbuFpEJjzcqFWMwY8raGwJPN53gYx/F2rwDU6yiuQisNfE8RX7Yi+YDB5l1uEC+cxcSjefMzGVA+9gjt1rU8M2+o21lKuoC4ViybVuJ/NbIjUOd2TwWDEDP4DOKANuiiigAooooA5GHw1fxRXsO2yKTiLd+8b/SGSUuTINn/AC0VirdcYH3geOg0eyk0/TI7aXYGVnYJGcpGGcsEXgcKCFHA4HQdKvUUAFFFFABWBc6Ndt4nh1OJbZ0WRWLu5WQJsZGjGFPy5bf15YYx3G/RQBjeH9Mn05bx7i2s7Z7ibzBHZuTGBjA42LzgDJ5yfQYA2aKKACiiigDN1vT5NRso4kjgmCTLI0FwcRygfwscH69DyBWRH4e1Qz+HZZby2I0sBXjZWfdiN4y6tkfMwZeCPl5wTznqaKACiiigAqK6iae0mhSQxPIjKsi9VJGMj6VLRQBxv/CK3r6U9jJb2S25uUn+yRXDpG4EWwoSEyo3ASZAOSSPc9XYwy29hbQzyCWaOJUeQLtDMAASB2ye1T0UAFFFFABXLvoF+NUvriIWm25iuIzM0jB38wKU3ADnYU2j5ujEggjB6iigDF8MaRPoumPbT+UCZS6iIggAgD+FEXOQeiqPqck7VFFABRRRQBk6/p9xqFrbLawWss0N3DOGuHKbAkisSpCsckAjt161Wi0a8PiFNTlW2UlxI7o5aQDydhgHyjMe7585HP8AD3rfooAKKKKACiiigAooooAKKKKAP//Z", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer: step: 1 , size: 64\n", "\t | u_h - u | = 1.92874e-05 , eoc = 3.01\n", "\t | grad(uh - u) | = 9.99164e-04 , eoc = 2.0\n" ] }, { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrASEDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3yV2jhd1ieVlGRGhGW9hkgfmRWTL4ltItOs73ybh47m2N3hVXMcIClnbLdt65AyeeAa1pY1mieNiwVwVJRip/AjkfUVmJ4a0pLOC0EMxggXZGr3MrYTABTJbJQhVyp+U45FAFm/1nS9KKDUdSs7Pf9z7ROse76biM1S/4TDwx/wBDHpH/AIHR/wDxVYXxNj83SdKT/qIj/wBETV5/9g9q8jHZp9Vq+z5b6X3Pay/KoYul7SU7a22/4J67/wAJh4Y/6GPSP/A6L/4qj/hMPDH/AEMekf8AgdF/8VXkBsfamS2P7p+P4TXGs/X8n4/8A7v9XqX/AD9/D/gnsX/CY+GP+hj0j/wOi/8AiqP+Ex8Mf9DHpH/gdF/8VXjSWP7pOP4RQbH2p/28v5Px/wCAUuHaT/5e/h/wT2X/AITHwx/0Mmj/APgdF/8AFUn/AAmXhf8A6GTR/wDwOi/+Krxc2PtUEdj/AKzj+M1SzxP7P4lf6uUr/wAX8P8Agnt//CZeF/8AoZNH/wDA6L/4qj/hMvC3/Qy6P/4HRf8AxVeINYe1RNY+1Us7X8v4lrhmk/8Al9+H/BPc/wDhM/C3/Qy6P/4HRf8AxVH/AAmfhb/oZdG/8D4v/iq8Eax/0jp/D/WmtYf7NX/bC/l/EqPC9J/8vvw/4J77/wAJp4V/6GbRv/A+L/4qk/4TXwp/0M2jf+B8X/xVfPzaeP7tVrjTx5bfLVrNk/s/iW+FKdr+2/D/AIJ9Ff8ACa+FP+hm0b/wPi/+Ko/4TXwp/wBDPov/AIHxf/FV85Npw7A1C2nHtVrNIvoP/VKPSt+H/BPpP/hNvCn/AEM+i/8AgfF/8VR/wm3hP/oZ9F/8D4v/AIqvmdtPYdqgSybyx8tWsxT6EPhLWyqfh/wT6f8A+E28J/8AQz6L/wCB8X/xVH/CbeE/+hn0X/wPi/8Aiq+XjaMP4ajNqfSmseuxL4Ukv+Xn4f8ABPqT/hN/Cf8A0NGi/wDgwi/+Ko/4Tfwn/wBDRov/AIMIv/iq+VjbfvOnag23tVfXV2Mv9WJfz/h/wT6p/wCE38J/9DRov/gwi/8AiqP+E38Jf9DRon/gwi/+Kr5UNt7VG9t04701jF2Ilw1Jfb/D/gn1f/wm/hL/AKGjRP8AwYRf/FUf8Jv4S/6GjRP/AAYRf/FV8o/Zvag23tT+uLsL/VuX8/4f8E+rv+E38Jf9DRon/gwi/wDiq2LS8tdQtY7qyuYbm3kGUlhcOjDpwRwa+OPs3tX1D8Mxj4c6KPSFv/Q2rajX9o2jzcxyt4KMZN3ubupa3pOjeV/amqWVj5ufL+1XCRb8YzjcRnGR+YqivjXwo7BV8TaMzE4AF/EST/31XmH7QC7p/Df+7dfzhryHT48alan/AKbJ/MUp1+WfLYrDZX7fD+25rb9Ox9eXWqwW9tb3KK1zbzzRwiWBlZV3sEVjkjIyQOMn2qleeKLGyv7m0kjuGa2Qs8ioNm4Jv8sEkZcqQcen0NaGoabbanAkN0JSiSLKBHM8Z3KcqcqQTggHHTiq8nh7TJmZpYHkLJsbfM7bvk2bjk8vt+Xf97Heug8gqx+KrKS+gszb3STyMUkDKuIDvZBvYNgbmUgYJz+Iq5ousR63Ym5S1urVlfY8F2gSWM4DYZcnHDA/jTF8PaYjRsIH3RncCZnJY7i4LfN85DEsC2cE5GKt2Wn2+npItur/ALx/MdpJWkZmwBksxJPCgde1AFmiiigAooooA5Tx1H5trpCf9P8A/wC0Jq5k2HtXY+KYxI2jqf8An+P/AKImql9jFfA8TVnDGpf3V+bPZwGKVKly36nNGw9qimsf3MnH8Jrqvsa+lRT2Q+zycfwn+VfPrFO53fX13OYjsf3Kcfwig2PtXSw2Q+zx8fwj+VKbH2p/W9TSOO8zljY+1VYbH/W8f8tDXYGx9qqQWOfO4/5atWkcXozRY7Xc5xrD2qJrH2rq2sfaomsfarWLNo47zOPax/0vp/B/Wlaw9q6RrH/TcY/5Z/1pzWHtW31vY1hjvM5NrD2qtdWOIG49P512DWPtVS8scWz8en861hi9UbfXtNzmWsfaomsfaurax9qiax9q0jizojjvM5RrH2qtFY5hXj1rr2sfaqsFjm3Xj1/nW0cXoarHa7nNNYe1RNYf7NdW1h7VE1h7VosWbxx3mci1gPPxt/h/rSNpy/3a6Y2P+lYx/B/WnGx9q1+t+ZccYjlDpw/u1BLpwGzj+IV15sPaq89j9zj+MVccXruW8XFo5n+zfakOnfWuq+w+1IbD/Zp/XPMv6zA5U6ca+hvhyuzwBpK+iOP/AB9q8h+we1exeAl2+CtPX08wf+RGr2Mpr+0nJeR8txTVjOlTt3f5Hnnx5XdceHP9y6/nDXlFlHi+tzj/AJar/OvX/jfH5l34eHpHdfzhry63t9lzExOArgkn61tiZ2r29DqyaipZVzf4v1PrCiqN1qsFvbW9yitc2880cIlgZWVd7BFY5IyMkDjJ9qgl12OHVZtPayvA6RGVJSirFNgAlUZmALc+2MGvWPgTVorn4PGGnzy2cflTp9r4jZjHjJdkABDHdll6ruGCpJAINaum366lZi5WGWH95JG0cu3crI7IwO0kdVPQ0AW6KKKACimSyxwRPLNIscaKWd3OAoHUknoKrvqunRwwTPf2qxTnbC5mULIfRTnn8KAMnxUu46ONzL/px5U4P+omqj5P/Tef/vur3is4Ojn/AKfj/wCiJqob6+Qz3DOriVLyX6nHXxPsp8o7yP8Ap4n/AO+6jng/0eT9/P8AdP8AH7U/fUc7/wCjyf7p/lXi/UXcyWO8x0EH+jx/v5/uj+P2qTyP+nif/vuooH/0eP8A3R/KpN9ZywTuWsd5jvs//Teb/vuoLW2B8799N/rW/jqbfUNq/wDrv+urVk8E7M2jjvMn+yKf+W0//fdJ9iQ/8t5/++6fvpQ9YSwkkbRx77lI6epv8efN/qv7/vUp0xT/AMtpv++6cH/4mH/bL+tWN9ZzoVFszohj2UW0sf8APab/AL6qpfaWBav+9m7fxe9bYeoL1gbR/wAP51EY1VJHRHHma2lf9NZv++qibSv+ms3/AH1W/wDKe1IVU1HtKqOmOOic42l/9NZv++qqW2l5tk/ey9/4veusMKmqtnbg2icev86tYiajqbxxke5z7aWP+ekv/fVRtpf/AE0l/wC+q6g2o9Kja09qaxTOiOLXc5E6Z/pmPMl/1f8Ae96edM/6aS/99V0BtP8AT8Y/5Zf1qQ2ftWrxb01NY4pHLnS/+mkv/fVV7jTMeV+8l/1g/irrTZ+1Vbq0/wBVx/y0FXDFu+5f1pWMD+zP+mkv/fVN/sz/AKaS/wDfVdN9j9qQ2ntT+tvuafWkcydM/wBuX/vqvR/BK7PCVmvPytKOf+ur1zv2P2rpvCI2+G4F9JZx/wCRXr6bhqs6lWovJHk5tW9pCK8zh/jAm+/0Ef8ATK6/nDXnSQEOp9CK9N+KieZqWhjHSG6/nDXCNBhScdBXZj6lsZb0PoMnq8uV8v8Ai/U981DTbbU4EhuhKUSRZQI5njO5TlTlSCcEA46cVEdGs2lWVvtLOsflqzXUp2jGMj5uGwPvDnk88mtCivpD4Mz4tEsIfI2Ry4gdpEVp5GXeWLFmBbDNlicnJzUun6bbaXbtBaiURtI0p8yZ5DuY5Y5ck8kk/iat0UAFFFFAFTUrL+0NPltfM8svgq+3IDAgjI7jIGR3rBuvCVxd2kcL6kilJZ5CyQMoPmtvYYEnOG3cNlSCAytjNdTRQBy/jBGEujS+c4UXpXysLtJ8iXnpnP449qz/ADK0PGxcQaR5e3d9v/i6f6iauf33X/TH9a5K+F9rLmPlM6xHs8Slfov1NDzKZNJ+4k/3T/Kqe+6/6Y/rTJnuvIk/1P3T6+lc/wDZ67HlLGeZfgk/cR/7o/lUu+suF7ryI/8AU/dHr6VJvuv+mP61nLL/ACLWM13NHfUNq/8Arv8Arq1Vd91/0x/WorZ7r97/AKn/AFh9axll5tHG6bmuHpd9Zwe6/wCmP607fdf9Mf1rnll5tHG+ZZ3/AOn/APbL+tWQ9Y++6+3f8sf9V7+tWA93/wBMf1rCeXnRHG+Zoh6hvH/0R/w/nVYPd/8ATH9ahu3uvsr/AOp7evrXM8v1N4Y3Xc2A9OD1nB7v/pj+tLvu/wDpj+tck8vOiON8zRD1DZP/AKIn4/zqqHu/+mP61FZPd/ZU/wBT39fWsJZfodEcb5mxvpd9Z2+7/wCmP60u+7/6YfrXNLLzeON8yxkHUOn/ACy/rVj5fSsnfd/b/wDlj/qvf1qxvu/+mH61lPAG0ca+5e2rVa6jX9z/ANdVqPfd/wDTD9ahuXu/3P8Aqf8AWr61ksC0zZY59zR8paTyV9qrb7v/AKYfrS77v/ph+tR9Sl3LWPfcseStXvCwxoKD0nuP/Rz1k77v/ph+ta3hXJ8Px7sbvPnzjp/rnr6jhag6dWo32X5iliHV0OV+JCb9V0Uf9MLn/wBChrjbiEJbSseAEJyfpXeePE36xo4/6d7n/wBChrlryDFlOfSNv5VeaVLZjb/D+h9Pl1XlwPL6nqN9qgt7O2u7VYbqCa4hiLrNgBZJFQMuAQ2Cw44+tVJvEcY8SnQ7dIJLkW7ykyT7PnGwhMYJ5V9xbt7841Lyws9QjWO9tILlEYOqzRhwGHQgHv70os7UIiC2hCRxmJF8sYVDjKj0XgcdOBX2R8qZOn67dXo0xpLGGNL+KSRGW4LYCnKH7gyGUg54IzjBq5o+oXGowXElxbRQ+XcPCpimMivs4Y5KqR8wZcY/hz3qw+nWMk8E72Vu01uMQyGJS0Q9FOMj8KmiijgjEcMaRxr0VFAA/AUAPooooAKKKKAOW8bsFt9IJ/5//wD2hNXPeaPWtr4hv5en6S3/AFEP/aE1cb9q969bBYX2tLm8z8/4n5vrqt/KvzZs+YPUU2aQeRJz/Cf5VlC696bLdfuX5/hNdLwB86pTubEL/uI/90fyqTfWNBdHyU5/hFTC7PrWUsAxupNM1N9Q2zf63/roaqi7pttdD97n/noawlgX2LjXkkzU30oeqguFPeniVT3rnngvIuOKY8N/p3/bP+tWA9UA/wDpvX/ln/WrAeuaWDOhYrYsB6iu2/0V/wAP50geortv9Ff8P51zSwZvTxWqLwenB6rB6dvrmng/I3jiywHqKzf/AEVPx/nTQ9RWb/6Kn4/zrmlgvI6YYvQvh6UNVcPS765pYI3jixd3+n/9sv61Y31Q3/6d/wBsv61Y31hLBG6xZY31Dct/qf8ArqtJvqG5f/U/9dVrJ4LXY1jiy/vo31Bvo31k8F5FLF+ZPvrW8K/8gCP/AK7z/wDo56wt9bnhP/kXov8ArtP/AOjnr1sooeynJ+R6WX1vaSkjJ8YJv1vShjpbXP8A6FDXP6hBjTLo46Qv/I11PiRPM17TB/063P8A6FDWRq0O3Rr5ueLeQ9P9k185nNS2a2/w/ofXYWry4fl9TvqKzr7U2trO2uobfzIpbiGFxKWiZRJIqA7SucgsODj61m6l4utrDWptMWOOWWK1eZszBTvGzCbcZwQ4JboB684/QDxjo6K5dvGlrDe2VpdJDDLNLMk+bgYhEfmjcMgFwTCw6DHfnAOl4f1yPX7CW6iREVJ3i2rKJOnTJHQ4IJHOPU9aANaiiigAoqlq8l1FpVw9mG88L8pRdzAZ5KrzkgZIGDkjoa5u41LX/sNqbcXe5JZxIzWR8yba48lCNuEDoeXwApHODkUAVvilKIdF0p2OB/aI/wDRM1ecjUYv7/6V3vxbklTR9KwieT9vGX3ndu8mXjGOnvn8K8v8/wB6+rySnzYZvzf6Hx+fUlPFJvsv1NT+0ov7/wChpJNRi8p/n/hPas4XGe9I9x+7bnsa9Z0HY8RYeN9jUh1GIRJ8/wDCOxqYalF/f/Q1kRT/ALteewqUXHvU+wZEsPG+xqjUov8Anp+hpLfUoh5nz/xnsazxce9Jbz/f/wB81nKg7mbw8bPQ211KL+/+hqRdTh/56foayVuPepVuPes5UDCVCPY0V1OL7ZnzP+WfofWra6pF/wA9P0NYaz/6V1/g/rVpbj3rnlhUyJ0UjXXVIO7/AKGkudRtzbP8/p2PrWetx70Tz5t259P51zTwSsZxjaSNkX9uf+Wn6Gni+g/56foazluPepluD61zTwIvaNFwXsH9/wDQ1FaXsH2ZPn9ex9aRbn1pLOZTbJ+P86454N9jWNe0S19ug/56foaX7dB/z0/Q00SA96cGrnlhDSOJRF9tg+3f6z/ln6H1qf7dB/z0/Q1Bu/03/tn/AFqxurnlhTf6ytA+3Qf89P0NQ3F9B+5/ef8ALQdjU+6obhv9V/10FZPCo1hiVcm+3Qf89P0NH26D/np+hpd1G6s3hV2H9aQfboP7/wChrp/CDBvDcDDoZZyP+/r1zO6um8I/8i3B/wBdZ/8A0a9XTo+zdz6Lh+rz1JryItbAPiDTc/8APrc/+hw1n6yF/sPUP+vaT/0E1f1048Qab/163P8A6HDWbrDf8SS//wCvaT/0E18LnNHmzXm/w/ofVqvyrlOsvLCz1CNY720guURg6rNGHAYdCAe/vQbG0MSxG1gMaRGFU8sYWM4BQD+6cDjpwKsUV9+SQR2drDHDHFbQokBLRKsYAjJBBKjscMRx6n1qSOKOLd5capuYs20YyT1J96fRQAUUUUAFFFFAHnXxkbZ4Z0w/9RJf/RM1ePC4z3r1v43Ns8J6af8AqJr/AOiZq8P8/wB6+64cp82Db/vP8kfP5pS5q1/I1vP96HnzG3PY1lifNDT/ACn6V7ro6HmewNaOf92vPYVKJ+etY0c/yjntUouPeo9joRKga4uPelhn+/z/ABGssXHvSwz/AHuf4jWcqOqM3Q0Ntbj3qVbj3rGWf3qVZ/eodExlQNVZ/wDSev8AB/WrSz+9YKz/AL/r/D/WrK3HvWPsTGdA2ln96Waf9w3Pp/Ospbj3p0s/7luf85rGdHRmPsNUbqz+9Srce9Yy3HvUyz+9ZyonNKgbK3HvRaT/AOjpz6/zrLWf3pbWf9wvP+c1hKjqZOh7rN1bj3qVbg+tYyz+9SrP71jLDJ9DB0WjTW5/00f9c/61aFwDWCJ/9L6/wf1qyLj3rmeETFKElY2BMpqOeQHyuf8AloKzxP70yaf/AFXP/LQVjPBCjzJm1vHrS7veswXGe9KLgjvWbwTFzSNLdXWeEP8AkWrf/rrP/wCjXrg/tHvXdeDDnwran1eb/wBGvXDjKHs0mfWcJybrVL9l+ZB4gONf0z/r1uf/AEKGsvV3xol+f+neTp/umtLxIca9pn/Xrc/+hQ1kaq3/ABJ73/r3k/8AQTXxGYUObHc3ofS4ity1+X0Ol1nWl07T7S9WeCCOW6hiYXaFCyvIqkAEqVYAluQenSq134gaPWYraC9sWtZbZpRIQG8vC7lYkScqevKquCPnyQD0VFfUHpnK2Wvapew6XJE1o7XEs4miWEkmONmBdWEhAziNf4hl+pFX/DOrz6xZSzTtbuVZAHt1IX5o1YryTypYg/yHStuigAooooAKKKKAPLvjuceDNNz/ANBRP/RM1eA7lr3r4+tt8Eacf+oon/omavnkTV+icK01LAttfaf5I87FwvO5dDL6frSll2n6VS873pTN8p57V9HKjG2xyezZeRl2jjt61KGX0/Ws5JflH0qQS+9R7GNtiJU2aAKen60sRT5uP4vWqQl96dHL1+tZyoxutDN03Y0lKen61KpT0/Ws9ZvepVlqHRj2MZU2XV2ed0/h9asLs9P1rMWX971/hqwstY+xWuhjOmzQXy/T9adIE8o8frVJZfenvL+6NYzoqz0MfZu5pqI/T9alUR+n61QWWpVl96zlRXY55U2aCrH/AHf1NLbLH5K/L+pqms1Ot5f3S81jKir7GMoOxpqsX939TUirF/d/U1RWX3qVZves3RXYwlBlkJF9q+7/AAep9asBIv7v6ms5Zf8ASf8AgH9asCWsvZLsZzgy4Ei/u/qaSVIv3fy/xjuariX3pJZf9Xz/ABis5UlbYzUJX3L4SL+7+pp2yL+7+pqp5tO873peyXYy5JFoJF/d/U16f4HwPCFljpul/wDRr15OJa9X8CnPg6xPvL/6MavFzmCjCPqfTcMxaq1L9kV/E5xrul/9e1z/AOhQ1jao3/Epvf8Arg//AKCa1vFZxrml/wDXtc/+hQ1h6m3/ABKbzgn9w/H/AAE18jWoc1XmOjMK3LjuX0PSaKxtZ1h9M0+0u3eC08y6hieO7xkq8iqQCHwGAJbPI4rO1PxS9teXSafc2F4sVkbnyVOGTIBQltx3AjLYC8LznkZ7z6k6qiuNTxVqe+w/c2dyJyyNHbEFpWWR1Zky+dgChgdrZGc7a1vC2sXOtaY9xdC381ZAubbJTBRXxkk5I3YPuDQBuUUUUAVNSvf7P0+W6EfmMmAqbtoLEgDJ7DJGT2rGvPF8NhY2ctxFEk9xeG1MbThVXbN5TuGIGQCQQMZOR05I6GWKOeJ4pY1kjdSro4yGB6gjuKhi06yggMMNnbxxF1cokShdwxg4A6jauD7D0oA8u/aBlUeDdNhw+46kjZ2HbjypR97GM89M5r53zX0V+0F/yImnY/6Ckf8A6Jmr5xyfav0ThWry4Fr+8/yRyV17xJmgng0zJ9qQlsdq+mdZW2MbEgJAFPDkVCCcdqXJ9qcaia2E0TiWnRy9ee9V8n2pF3c9OtDabWhLgi+stSrL71mhnHpT1kf2ocfIzdI0ll/efhVhZveslZG39R0qZZH/ANmsOXfQxlSNVZae0v7s1mLI/wDs1IZH2H7tZVI+69DF0tTXWX3qZZqyVlk/2amWST/ZrOUfI55UTVWb3p0Ev7pazVlk/wBmnQySeWOVrCUddjGVHQ2Fl96lWaspZZf9j9alWWX/AGP1qHDyMJUTRWX/AEj/AID/AFqcTe9ZAkl87+D7vv61MJJfVP1rHk8jOVE1RNSSS/6vn+MVnCWb/Y/Wh5Zfk+594etRKGhn7HU1xL707zveswSzf7H60olm/wBj9aXIZ+xNMS17H4COfBWnn/rp/wCjGrwoSzeqfrXufw+JPgbTCeuJP/RjV8/n8bU4ep72QQ5ak/QreLjjW9K/69rn/wBChrB1Fv8AiWXf/XF//QTW140cJrWkk/8APvc/+hQ1zl/MDp1yPWJ/5GvCp4fnhzHkZzVcc15f8P6Hq1FFFcp98FFFFABRRRQAUUUUAeT/ALQf/Iiad/2FI/8A0TLXzhmvo79oT/kQ9P8A+wrH/wCipa+b8191w3PlwbX95/kjlrfEOzQTwabmgng19C6mhlYeDwKXNMB4FLmtIVdAsPzSKev1pM0Kev1rVVNUKw/NLTM0ua6YzTJsOBw/4VKslQ5+b8KWnGzuJq5bWSnl/lNUwxFP35U1FSPuszcDQWSplkqgslTK9RKBzygX1kp8L/uxVJZKkik+QVzyh7xjKGhoLJUivVJZKkV6hwMJQLayfvv+A1OJKz1k/e/8BqYPWKgZSgXRJQ7/AHP94VWD0PJ9z/eFROGhnyal4PSh6qiSlD+9LkM+Qth695+HvPgTTP8Adf8A9GNXz6Hr6A+Hf/IhaV/uP/6MavmuJI2pQ9Wezk0bTl6GV8QZPL1bRz6wXP8A6FDXI3VxmznABP7tuPwrpPidJ5ep6Ic9Ybr/ANChriZLjdE4z1UiuXAU74O/qeJnFLmzTm/w/oeyaze3UGn2lzAl7DKbqEPDFAJmKGRQ4YKHwNm45BHQc9qz9R1LUmv7gaa2obFszJtksiE3EAjyyYwS4GSQxPJC4Jzt6mivmz7s5ODUNTa7sYpJtUAJmaV2sPkaENIIycR5EhG04yBgZI5AbW8O3F1daaZrprolpGMYu4PKlVOMBxtUZ78DHOMnGTrUUAFFFFABRUVzcw2dtJcTvsijGWOCfyA5J9h1qm2uWCWsVyZJfJlcoHFvIQrBtpD/AC/JhuDux0PpQB5t+0KceA9P/wCwrH/6Kmr5t3fWvpL9oX/kQtP/AOwrH/6Kmr5tzX2GQyawr16v9Dnq/ELu9jQW46GkoPSvb5nbczFDcdDS7vY00dKXNaQk7biHbvY0BuvBpKB3+taqTutRD93saXd7Gm0tdMZPuIXd83Q9Kdu9jTP4vwp1b05S7iY7d7GkZvlPBpaD92uh3cXqSPEhHY1KsvsagpwOKpxfclpMtLN7H8qkil+UcN+VV1anxNhRWEoPm3MZRRbWb2b8qlWb/Zb8qrK1SK1Q4PuYSiicTfvPut930qYTf7LflVVW/efhUwasVB66mcoonE/+y35UrT/d+VvvDtUQalZvu/7wqJwdtzPlVywJ/wDYb8qUT/7L/lUQalDUuR9zPlRL5/8Asv8AlX0Z8ODn4f6SfWN//Q2r5wDV9H/Df/kn2kf9c2/9Davl+J01Sp+r/I9PK1acjmPi3J5eoaEfWK6/nDXnwuNxA9eK7f40SeXe+Hz6x3X84a8ziuCZkA5JYYzTyqnfLr/4v1ODMKXNjub0PqaisbWbu9tdPtJ0E8c/2qFZIrSI3AZDIocH5Cduzcc4Xp1qKa/nOsuiSagtpLYGRcWTYjfIxg+WTuxn5Tnp93tXxh9Sb1Fcbb6j4lJ0jfaXjwmaYXMoijy6/vdgYMEKgbYznau7I6dDp+EbrVbrS5W1fzvPEoCmWLyyV8tCeCiZ+Yt/D7ZOMkA36KKKAK99Zx39nJbSsyq4HzIfmUg5BGe4IB/Csz/hHE8qCL+077y4blrrYfKxI7NuO4bORuJYDsTx0XG3RQB5D+0NEv8AwhOnTZfcNSRcbztx5Up+7nGffGa+cM19JftDf8iFp/8A2FY//RUtfNma+ryWVsM/V/oYVNxaD0pKO1exzaGYo6UuaaOlLWkJ6CHUDv8AWkoHetlLVAPpabS10xkSL/F+FOpmfm/CnVvTkJjqUnim0p6V083uskfRSUtdSdxDlbFSxt8oqClQ4UVLXvIlq6LitUitVVWqVWqGjGUSwrfP+FShqqhvn/CpQ1YpbmUolkNSs33f96oQ1Kzfd+oqJrQy5dSyGp26oN1ODUrEcpMGr6T+G3/JPdH/AOubf+htXzOGr6X+Gv8AyTzR/wDrk3/obV8pxUrUqfq/yPRy5WlI4b46yeXc+HTnql1/OGvKrS43XsAz1kUfrXpf7QUnlzeGz6rdf+0a8g06fOp2g9Zk/mK7Mlp3ym/lL9QxNLmr83ofZtFFFfnp6wUUUUAFFFFABRRRQB5J+0P/AMiDp/8A2FY//RU1fNWBX0r+0P8A8iDp/wD2FY//AEVNXzVX02T2+rv1/wAjGpuGBRgYpaD0r1rKxmAAxS4FIOlLWkUrbALgUACigVqkrrQQuBS4FFLXTFLsSGBupwApvf8ACnVvTUewmKAKCBiig9K6eWPK9BDsClwKBS11wjHsSG0elIqjHSnUi9BWnJHmWgh6qKlCr6VEDUqmhwj2IkPVV8zp2qUKvpUQ+/8AhUoNYqEddDKQ8KvpTiqkDjvTRTifu/WonCNtjLUeEX0pdi+lFKDS5I9iHcUIp7V9O/DTj4daN/1yb/0Nq+Yq+nfhp/yTrRv+uTf+htXyfFcUqVO3d/kd2A+KR5v+0UcSeGfpd/8AtGvGdMY/2rZ45PnpjJ/2hXsn7RhxJ4Z+l3/7RrxjSm/4m9l/13j/APQhXXktblynl/xfqdNRfvD7D1ldQk0+0eC1na7S6hd0tLkKAgkUvksUDLsDDBHfp3rM1O18QXNxdS2QuoJJI2MTNOuxEMGBHsDEeYJfm3Yxj+LHFdZRX56dZy9rpeqFdM8xr1Wt5ZpHaa9bLoC3lpIFYhiSyknBGEI7gVe8OW+rW0F4ursGla5LxsJzICpRc7cgbV3bsL2raooAKKKKACiqWrx3UulXEdmWE5X5QjbWIzyFbjBIyAcjBPUVzs1pr0mnRRxx3yOhuVjzdruVmYGB3YP86IpKsCWJI6NwaAOT/aH/AORB0/8A7Csf/oqavmqvr/4heCpPHWnWOmPei2tIroTysse58iN1GCTjGWAxj8a4T/hnLTf+hhu/+/C/417WX46lQpck97mc4ts+e6O1fQn/AAzlpv8A0MN3/wB+F/xo/wCGctN/6GG7/wC/C/413f2th/P7iPZs+fB0FLX0F/wznpv/AEMN3/34X/GmT/s62SW8jQa9cyShCURoVAZscAnPHNXHOMMu/wBwezZ4DQK+gv8AhnTTf+hhu/8AwHX/ABo/4Z000f8AMw3f/gOv+NaLOsLfr9wvZyPn4U4V9Af8M66d/wBDDd/+A6/40f8ADOunf9DDdf8AgOv+NbRz3CLq/uF7KR8/j71Or37/AIZ107Of+Ehuv/Adf8aZH+zxZl5RJrtyqh8RkQqdy7RyeeOcj8K1hxBg1u39weykeC0HpXv3/DO+nf8AQwXX/gOv+NH/AAzvp+P+Rguv/Adf8a2/1jwVrXf3C9jI8DpRXvn/AAzzp/8A0MF1/wCA6/40v/DPOn/9DBdf+A6/410R4nwC6v7ifYyPBBSL92veD+z5aC4RRrtwYijFn8lchsjAxnuC35e9SD9nrTwP+Rguv/Adf8a0/wBacvve7+4PYTPBRUiGvdv+GerD/oYLr/wHX/GnD9nywH/MwXP/AIDr/jQ+Ksv7v7iXh5s8LB+f8KlWvcP+GfrHOf8AhILn/wAB1/xpsvwDt0QGHXJ3beoIMCj5Sw3Hr2GT+FZrifL+7+4h4aozxQGlPb617f8A8KDs/wDoP3H/AIDr/jQfgHZnH/E/uOP+ndf8amXE2Aa3f3EfVKh4mDing5r2v/hQ1p/0H5//AAHX/wCKo/4UNaf9B+4/8B1/+Ko/1mwHd/cT9TqnigNfT/w0/wCSdaN/1yb/ANDauM/4UPaf9B+f/wABl/8Aiq9F8KaPJ4f8N2ulSSCQ2xkRXAxuXexUkeu0jNeBn2a4fHU4Ro3un1R04WhOk25HkP7R5w/hj6Xf/tGvFtJb/ic2P/XxH/6EK+qviD8OLb4gHTjcahNZ/YfN2+XGG3b9mc59Ng/OuOt/2ebC2uYp4/EFyXicOu63UjIORnmufB5nTo4X2Mt9fxOiUG5XPZ6KxtZi1A6faGCOe5u47qF3NpIIQUEil8hnAKlNw2knrWZqNrr1zdXMtqt3C0iEx5uVCrGYMeVtDYEnm87wMY/i7V4BqdZRXIRWGvieIr9sRfMBg8y63CBfOYuJRvPmZjKgfewR261qeGbfUbaylXUBcKxZNq3E/mtkRqHO7J4LBiBn8BnFAG3RRRQAUUUUAcjD4av4or2HbZFJxFu/eN/pDJKXJkGz/lorFW64wPvA8dBo9lJp+mR20uwMrOwSM5SMM5YIvA4UEKOBwOg6VeooAKKKKACsC50a7bxPDqcS2zosisXdysgTYyNGMKfly2/rywxjuN+igDG8P6ZPpy3j3FtZ2z3E3mCOzcmMDGBxsXnAGTzk+gwBs0UUAFFFFAGbrenyajZRxJHBMEmWRoLg4jlA/hY4P16HkCsiPw9qhn8Oyy3lsRpYCvGys+7Ebxl1bI+Zgy8EfLzgnnPU0UAFFFFABUV1E09pNCkhieRGVZF6qSMZH0qWigDjf+EVvX0p7GS3sltzcpP9kiuHSNwIthQkJlRuAkyAckke56uxhlt7C2hnkEs0cSo8gXaGYAAkDtk9qnooAKKKKACuXfQL8apfXEQtNtzFcRmZpGDv5gUpuAHOwptHzdGJBBGD1FFAGL4Y0ifRdMe2n8oEyl1ERBABAH8KIucg9FUfU5J2qKKACiiigDJ1/T7jULW2W1gtZZobuGcNcOU2BJFYlSFY5IBHbr1qtFo14fEKanKtspLiR3Ry0gHk7DAPlGY93z5yOf4e9b9FABRRRQAUUUUAFFFFABRRRQB//9k=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dt.value = 0.005\n", "\n", "errors = 0,0\n", "loops = 2\n", "for eocLoop in range(loops):\n", " u_h.interpolate(initial)\n", " evolve(scheme, u_h, u_h_n, endTime)\n", " errors_old = errors\n", " errors = [sqrt(e) for e in integrate([l2error,h1error])]\n", " if eocLoop == 0:\n", " eocs = ['-','-']\n", " else:\n", " eocs = [ round(numpy.log(e/e_old)/numpy.log(0.5),2) \\\n", " for e,e_old in zip(errors,errors_old) ]\n", " print('Forchheimer: step:', eocLoop, ', size:', gridView.size(0))\n", " print('\\t | u_h - u | =', '{:0.5e}'.format(errors[0]), ', eoc =', eocs[0])\n", " print('\\t | grad(uh - u) | =', '{:0.5e}'.format(errors[1]), ', eoc =', eocs[1])\n", " u_h.plot()\n", " if eocLoop < loops-1:\n", " gridView.hierarchicalGrid.globalRefine(1)\n", " dt.value /= 2" ] }, { "cell_type": "markdown", "id": "d2b700c0", "metadata": {}, "source": [ ".. index::\n", " pair: Solvers; Scipy\n", "\n", "## Using Scipy\n", "We implement a simple Newton Krylov solver using a linear solver from\n", "Scipy. We can use the `as_numpy` method to access the degrees of freedom as\n", "Numpy vector based on the `python buffer protocol`. So no data is copied\n", "and changes to the dofs made on the Python side are automatically carried\n", "over to the C++ side." ] }, { "cell_type": "markdown", "id": "8652a2e8", "metadata": {}, "source": [ "The most important step is accessing the data structures setup on the C++\n", "side in Python. In this case we would like to use the underlying dof vector from\n", "a discrete function as numpy arrays and system matrices assembled by the\n", "schemes and operators as scipy sparse matrices.\n", "In the [introduction](dune-fempy_nb.ipynb) we already discussed the `as_numpy` method.\n", "So" ] }, { "cell_type": "code", "execution_count": 5, "id": "3627dc88", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:10:47.891611Z", "iopub.status.busy": "2024-03-01T11:10:47.891386Z", "iopub.status.idle": "2024-03-01T11:10:47.896409Z", "shell.execute_reply": "2024-03-01T11:10:47.895060Z" } }, "outputs": [], "source": [ "vecu_h = u_h.as_numpy" ] }, { "cell_type": "markdown", "id": "330b6fb7", "metadata": {}, "source": [ "provides access to the underlying dof vector without copying. So changes\n", "to the numpy array `vecu_h` carries over to the discrete function `u_h`.\n", "Just remember to make changes using `vecu_h[:]` to change the actual\n", "memory buffer.\n", "\n", "A `scheme` describing an operator `L`\n", "provides a method `linear` which returns an object that\n", "stores a sparse matrix structure. The object describes the operator\n", "linearized around zero.\n", "To linearize around a different value use the `jacobian` method on the `scheme` that linearized\n", "the the operator `L` around a given grid function `ubar` and fills the\n", "linear operator structure passed in as second argument. It is also\n", "possible to pass ``assemble=False`` to the ``linear`` method to avoid an\n", "the linearization around zero to reduce computational cost:" ] }, { "cell_type": "code", "execution_count": 6, "id": "19cdc408", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:10:47.901156Z", "iopub.status.busy": "2024-03-01T11:10:47.900932Z", "iopub.status.idle": "2024-03-01T11:11:08.292878Z", "shell.execute_reply": "2024-03-01T11:11:08.291337Z" } }, "outputs": [], "source": [ "linOp = scheme.linear() # linearized around 0\n", "linOp = scheme.linear(assemble=False) # empty (non valid) linear operator\n", "scheme.jacobian(space.zero, linOp) # linearized around zero" ] }, { "cell_type": "markdown", "id": "69595d7f", "metadata": {}, "source": [ "Here we linearize around zero. But that argument could be any grid\n", "function. A second version of this method will return an addition\n", "discrete function `rhs` which equals `-L[ubar]` such that `DL[ubar](u-ubar) - rhs` are the first\n", "terms in the Taylor expansion of the operator `L`:" ] }, { "cell_type": "code", "execution_count": 7, "id": "e1b73e4c", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:11:08.299680Z", "iopub.status.busy": "2024-03-01T11:11:08.299110Z", "iopub.status.idle": "2024-03-01T11:11:08.307812Z", "shell.execute_reply": "2024-03-01T11:11:08.306577Z" } }, "outputs": [], "source": [ "rhs = u_h.copy()\n", "scheme.jacobian(u_h, linOp, rhs)" ] }, { "cell_type": "markdown", "id": "83e2908e", "metadata": {}, "source": [ "One can now easily access the underlying sparse matrix by again using\n", "`as_numpy` (and again the underlying data buffers are not copied):" ] }, { "cell_type": "code", "execution_count": 8, "id": "10297b16", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:11:08.313583Z", "iopub.status.busy": "2024-03-01T11:11:08.313058Z", "iopub.status.idle": "2024-03-01T11:11:08.319591Z", "shell.execute_reply": "2024-03-01T11:11:08.318559Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "A = linOp.as_numpy\n", "print(type(A))\n", "# plt.spy(A, precision=1e-8, markersize=1)" ] }, { "cell_type": "markdown", "id": "cb9cb8c0", "metadata": {}, "source": [ "Now we have all the ingredients to write a simple Newton solver to solve\n", "our non-linear time dependent PDE." ] }, { "cell_type": "code", "execution_count": 9, "id": "35c541f0", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:11:08.325810Z", "iopub.status.busy": "2024-03-01T11:11:08.325281Z", "iopub.status.idle": "2024-03-01T11:11:50.744796Z", "shell.execute_reply": "2024-03-01T11:11:50.743228Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer(numpy) size: 64 L^2, H^1 error: 1.93091e-05, 9.99162e-04\n" ] } ], "source": [ "import numpy as np\n", "from scipy.sparse.linalg import spsolve as solver\n", "class Scheme:\n", " def __init__(self, scheme):\n", " self.model = scheme.model\n", " self.jacobian = scheme.linear()\n", "\n", " def solve(self, target):\n", " # create a copy of target for the residual\n", " res = target.copy(name=\"residual\")\n", "\n", " # extract numpy vectors from target and res\n", " sol_coeff = target.as_numpy\n", " res_coeff = res.as_numpy\n", "\n", " n = 0\n", " while True:\n", " scheme(target, res)\n", " absF = numpy.sqrt( np.dot(res_coeff,res_coeff) )\n", " if absF < 1e-10:\n", " break\n", " scheme.jacobian(target,self.jacobian)\n", " sol_coeff -= solver(self.jacobian.as_numpy, res_coeff)\n", " n += 1\n", "\n", "scheme_cls = Scheme(scheme)\n", "\n", "u_h.interpolate(initial) # reset u_h to initial\n", "evolve(scheme_cls, u_h, u_h_n, endTime)\n", "error = u_h - exact_end\n", "print(\"Forchheimer(numpy) size: \", gridView.size(0), \"L^2, H^1 error:\",'{:0.5e}, {:0.5e}'.format(\n", " *[ sqrt(e) for e in integrate([error**2,inner(grad(error),grad(error))]) ]))" ] }, { "cell_type": "markdown", "id": "0176d448", "metadata": {}, "source": [ "We can also use a non linear solver from the Scipy package" ] }, { "cell_type": "code", "execution_count": 10, "id": "4aab17df", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:11:50.761028Z", "iopub.status.busy": "2024-03-01T11:11:50.760131Z", "iopub.status.idle": "2024-03-01T11:11:53.380487Z", "shell.execute_reply": "2024-03-01T11:11:53.378966Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer(scipy) size: 64 L^2, H^1 error: 1.93091e-05, 9.99162e-04\n" ] } ], "source": [ "from scipy.optimize import newton_krylov\n", "from scipy.sparse.linalg import LinearOperator\n", "from scipy.sparse.linalg import cg as solver\n", "\n", "class Scheme2:\n", " def __init__(self, scheme):\n", " self.scheme = scheme\n", " self.model = scheme.model\n", " self.res = u_h.copy(name=\"residual\")\n", "\n", " # non linear function\n", " def f(self, x_coeff):\n", " # the following converts a given numpy array\n", " # into a discrete function over the given space\n", " x = space.function(\"tmp\", dofVector=x_coeff)\n", " scheme(x, self.res)\n", " return self.res.as_numpy\n", "\n", " # class for the derivative DS of S\n", " class Df(LinearOperator):\n", " def __init__(self, x_coeff):\n", " self.shape = (x_coeff.shape[0], x_coeff.shape[0])\n", " self.dtype = x_coeff.dtype\n", " x = space.function(\"tmp\", dofVector=x_coeff)\n", " self.jacobian = scheme.linear()\n", " self.update(x_coeff,None)\n", " # reassemble the matrix DF(u) given a DoF vector for u\n", " def update(self, x_coeff, f):\n", " x = space.function(\"tmp\", dofVector=x_coeff)\n", " scheme.jacobian(x, self.jacobian)\n", " # compute DS(u)^{-1}x for a given DoF vector x\n", " def _matvec(self, x_coeff):\n", " return solver(self.jacobian.as_numpy, x_coeff, tol=1e-10, atol=1e-10)[0]\n", "\n", " def solve(self, target):\n", " sol_coeff = target.as_numpy\n", " # call the newton krylov solver from scipy\n", " sol_coeff[:] = newton_krylov(self.f, sol_coeff,\n", " verbose=0, f_tol=1e-8,\n", " inner_M=self.Df(sol_coeff))\n", "\n", "scheme2_cls = Scheme2(scheme)\n", "u_h.interpolate(initial)\n", "evolve(scheme2_cls, u_h, u_h_n, endTime)\n", "error = u_h - exact_end\n", "print(\"Forchheimer(scipy) size: \", gridView.size(0), \"L^2, H^1 error:\",'{:0.5e}, {:0.5e}'.format(\n", " *[ sqrt(e) for e in integrate([error**2,inner(grad(error),grad(error))]) ]))" ] }, { "cell_type": "markdown", "id": "982a4f8d", "metadata": {}, "source": [ ".. index::\n", " triple: Solvers; Dirichlet; Conditions\n", "\n", ".. index::\n", " pair: Boundary; Linear Solvers\n", "\n", "## Handling Dirichlet boundary conditions\n", "We look at a simple Poisson problem with Dirichlet BCs\n", "to show how to use external solvers like the cg method\n", "from Scipy in this case.\n", "We solve $-\\triangle u=10\\chi_\\omega$ where $\\chi_\\omega$ is a characteristic\n", "function with $\\omega=\\{x\\colon |x|^2<0.6\\}$. For the boundary\n", "we prescribe trivial Neuman at the top and bottom boundaries\n", "and Dirichlet values $u=-1$ and $u=1$ at the left and right\n", "boundaries, respectively.\n", "We will use the CG solver from `scipy.sparse.linalg`.\n", "\n", ".. tip:: Since we are not needing to invert the operator\n", "we will use the `dune.fem.operator.galerkin` class\n", "to setup the problem. This is similar to `dune.fem.scheme.galerkin`\n", "we have been using so far but can be used to model\n", "operators between different spaces.\n", "See [here](scheme_api.rst) for a summary of the concepts and API for\n", "operators and schemes." ] }, { "cell_type": "code", "execution_count": 11, "id": "132365ea", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:11:53.386299Z", "iopub.status.busy": "2024-03-01T11:11:53.385727Z", "iopub.status.idle": "2024-03-01T11:13:21.200112Z", "shell.execute_reply": "2024-03-01T11:13:21.197830Z" } }, "outputs": [], "source": [ "from dune.ufl import DirichletBC\n", "from dune.fem.operator import galerkin\n", "from scipy.sparse.linalg import cg as solver\n", "model = ( inner(grad(u), grad(v)) -\n", " conditional(dot(x,x)<0.6,10.,0.) * v ) * dx\n", "dbcs = [ DirichletBC(space,-1,1), DirichletBC(space, 1,2) ]\n", "op = galerkin([model, *dbcs], space)\n", "sol = space.interpolate(0, name=\"u_h\")\n", "rhs = sol.copy()\n", "lin = op.linear()" ] }, { "cell_type": "markdown", "id": "a8328a62", "metadata": {}, "source": [ "So far everything is as before. Dirichlet boundary conditions\n", "are handled in the matrix through changing all rows\n", "associated with boundary degrees of freedom to unit rows -\n", "associated columns are not changed so the matrix will not be symmetric anymore.\n", "For solving the system we need to modify the right hand side\n", "and the initial guess for the iterative solver to include\n", "the boundary values (to counter the missing symmetry).\n", "We can use the first of the three versions of the\n", "`setConstraints` methods on the scheme class discussed in the section on\n", "[more general boundary conditions](boundary_nb.ipynb#Accessing-the-Dirichlet-degrees-of-freedom)." ] }, { "cell_type": "code", "execution_count": 12, "id": "d9028d91", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:13:21.208074Z", "iopub.status.busy": "2024-03-01T11:13:21.207485Z", "iopub.status.idle": "2024-03-01T11:13:21.594702Z", "shell.execute_reply": "2024-03-01T11:13:21.593363Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "22.364801469753214\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "6.5503071305549865\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "5.502038979007685\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2.6517035748641966\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2.8478465378604483\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "1.3560658457805435\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "1.6365176519674627\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "1.199620085166859\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.8491708284280252\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.41168512215051617\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.3318991774651992\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.25180064032420496\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.236151715398552\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.22530579266807343\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.2214195904289019\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21926871065555859\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.2182915353763845\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21653684764907855\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.2167578210704593\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21581106387180027\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21559315250503042\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21544432425537105\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21544150394440711\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.2153861339325573\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538859378926486\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538505227526789\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.2153841894750399\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538241666417574\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.2153828225690404\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538287372050158\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538276135478884\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538277558134702\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538281353457003\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283870317615\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283905213898\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283539716058\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283650631818\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.215382836937197\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283753480741\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283758833887\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283758714316\t" ] }, { "name": "stdout", "output_type": "stream", "text": [ "0.21538283755554374\t" ] }, { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADrASIDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3yV2jhd1ieVlGRGhGW9hkgfmRWTL4ltItOs73ybh47m2N3hVXMcIClnbLdt65AyeeAa1pY1mieNiwVwVJRip/AjkfUVmJ4a0pLOC0EMxggXZGr3MrYTABTJbJQhVyp+U45FAE+pa5pOjBDqmqWViH4U3VwkW76biM1n/8Jz4R/wChq0T/AMGEX/xVcH+0P/yIWn/9hWP/ANFTV83LXsZflaxcOdytr2M5z5WfZv8AwnPhH/oadE/8GEX/AMVS/wDCceEv+hp0T/wYRf8AxVfGy1J/Ca92nwlCav7X8P8AgmTr+R9ijxv4SPTxRon/AIMIv/iqX/hNfCh/5mbRv/A+L/4qvka1+6tbFtRPhKEf+Xv4f8E5quPcPsn1H/wmfhY9PEuj/wDgdF/8VSjxj4YPTxJo5/7fov8A4qvnG37Vp2X8f++a46nDsYf8vPw/4Jw1M7lD7H4/8A98/wCEu8Mn/mYtJ/8AA2P/AOKpf+Es8Nn/AJmDSv8AwNj/AMa8dt+ta1v2rjqZMofb/D/gnn1OKZQ/5dfj/wAA9N/4Svw6Tj+39K/8DI/8aX/hKfD3/Qe0v/wMj/xrgoP+Pr/tn/Wr615Vaj7PqKHFMpf8uvx/4B13/CU+Hv8AoPaX/wCBkf8AjR/wlPh4DJ17S/8AwMj/AMa5daS6/wCPV/w/nXl1cU6fQ9Gjnkqn2Px/4B1X/CUeH/8AoO6Z/wCBcf8AjR/wk+gf9BzTP/AuP/GsJamWvKrZ5Kn9j8f+AejSxrn9k1/+Em0D/oOab/4Fx/40DxPoBGRrmmY/6+4/8azlqK0/48U/H+dedU4rlD/l1+P/AAD0aMfaM1v+Eo8Pjrrumf8AgXH/AI03/hKvDo/5j2l/+Bkf+NYdx1rKuKUOK5S/5dfj/wAA9Knlyn9r8DsP+Es8Ng4/4SDSs/8AX7H/AI0n/CXeGh/zMOk/+Bsf/wAVXmtx/wAfh/65/wBaoXHWuyPEMpf8u/x/4B3UsijP7f4f8E9Y/wCEw8Mjr4j0j/wNj/8AiqQ+MvC46+JNHH/b9F/8VXidx0rHvf4f94V1084c/sfj/wAA7ocLxkv4v4f8E+gz408Kjr4m0b/wPi/+KpP+E28Kf9DPov8A4Hxf/FV80XX3qzZeldcMe5fZHV4VjT/5e/h/wT6o/wCE38J/9DRov/gwi/8Aiq2LS8tdQtY7qyuYbm3kGUlhcOjDpwRwa+Mpe9fU3ww/5Jton/XA/wDobV10qvtOh4OYZesIlaV7m9qWt6To3lf2pqllY+bny/tVwkW/GM43EZxkfmKor418KOwVfE2jMxOABfxEk/8AfVeYftALun8N/wC7dfzhryHT48alan/psn8xUzr8s+Wxphsr9vh/bc1t+nY+vLrVYLe2t7lFa5t55o4RLAysq72CKxyRkZIHGT7VUvfElpYXE8c0NwVhJRpVVSpkEfm+WOc7tnPTHbOavahpttqcCQ3QlKJIsoEczxncpypypBOCAcdOKryeH9MmZmlgeQtH5bb5nYN8mzccnl9vG/72O9dB5BXj8TWsk6wm2ukYOsc24J+4ZpDGgbDc7mUgbd3qcCtusyPQNNjkikWB90TbgTM53HcXBfJ+chiSN2cE5GK06ACiiigAooooA8k/aH/5ELT/APsKx/8Aoqavm5a+kf2h/wDkQtP/AOwrH/6Kmr5uWvrsg/gv1/yMKu5KtSfwmo1qT+E19vhvhOZl21+6tbFtWPa/dWti2pVdjzcSa1v2rTsv4/8AfNZlv2rTsv4/9815Ffc8Svszbt+ta1v2rJt+ta1v2rya+x4OINGD/j6/7Z/1q+tUIP8Aj6/7Z/1q+tfLYvdio7omWkuv+PV/w/nSrSXX/Hq/4fzr5vFbnvYTdF1amWoVqZa+XxZ9DhiZaitP+PFPx/nUq1Faf8eKfj/OvncRsz6HC7or3HWsq4rVuOtZVxSon0eHMS4/4/D/ANc/61QuOtX7j/j8P/XP+tULjrXr0uh7WGMq46Vj3v8AD/vCti46Vj3v8P8AvCvUobns0tjKuvvVmy9K0rr71ZsvSvVpbGmK2KUvevqb4Yf8k20T/rgf/Q2r5Zl719TfDD/km2if9cD/AOhtXq4XqfC59tH1OF+PK7rjw5/uXX84a8oso8X1ucf8tV/nXr/xvj8y78PD0juv5w15db2+y5iYnAVwST9a5cTO1e3oe1k1FSyrm/xfqfWFFZeoa/Z6fax3ID3MTo8oa3KsBGgyz5JAIA9MnngGn3GrG31RLAWF1LJJE8qOhj2sExkcuCDllHIAyetesfAmjRXPw+L7OW4sYGtbqKS9ZhEshjzw5TIAc7h8ucpuwpBOAa6CgAooooAKKZLLHBE8s0ixxopZ3c4CgdSSegqu+q6dHDBM9/arFOdsLmZQsh9FOefwoA8v/aH/AORB0/8A7Csf/oqavm1R7mvpL9of/kQtP/7Csf8A6Kmr5uWvrcgSdF+v+RhV3JFHuakx8p5NMWpP4TX22HiuU5mXLVflX5j+da1un+0351lWv3VrYtqKsVY87EM07ePp87/nWlZRff8A3j/fPeqFv2rTsv4/9815FeKueJXbszWt4ef9ZJ/31WtbwdP3sv8A31Wdb9a1rftXk10rHg4iTLcFv/pX+ul/1f8Ae96vrbf9Npv++qrQf8fX/bP+tX1r5fF7sVGTuhFtv+m03/fdJc23+iv++m7fx+9WFpLr/j1f8P5183itz3sK3dEy23/Teb/vupltf+m8/wD33SrUy18ximz6DDMatr/03n/77qK0tf8AQk/fz9/4/erq1Faf8eKfj/OvnsRJ2Z9Dhd0ULi25/wBfN/33WVcW/wD02m/76rcuOtZVxU0ZM+iw6OfuIP8ATD+9l/1f973qhPD/ANNZP++q1Lj/AI/D/wBc/wCtULjrXr0pPQ9nDpGNPDx/rJP++qyLyL7vzv8AeHetu46Vj3v8P+8K9Sg3c9qlFWMe6j+b77fnWdKn+0351qXX3qzZelerSeheKirFGVevJ/Ovqf4Yf8k20P8A64H/ANDavlmXvX1N8MP+SbaJ/wBcD/6G1ephep8Nnyso+pzHxgTff6CP+mV1/OGvOkgIdT6EV6b8VE8zUtDGOkN1/OGuEaDCk46CvHx9S2Mt6Hs5PV5cr5f8X6nveoaXZ6pGsd5EZFXOMOy5BGCDtIyCOoPB70+PT7eOXzVVzJ5Kwb2kZm2DOBknOeeT1PGScCrNFfSHwZm/2Bpge2Zbcr9mSOONVlcLtjOUBUHDbTyM5wa0qKKACiiigCpqVl/aGny2vmeWXwVfbkBgQRkdxkDI71g3XhK4u7SOF9SRSks8hZIGUHzW3sMCTnDbuGypBAZWxmupooA8h/aHjb/hCdOk819o1NF8vA2k+VLz0zn8cc9K+cVr6R/aH/5EHT/+wrH/AOipq+bVz7V9bkD/AHL9f8jCruTLUn8JqJc+1Sc7T0r7bDy93Y5mXrX7q1sW1Y1ru2r0rWt/M/2aKstNjzsQjZt+1adl/H/vmsm383j7n61pWXnfPjZ98+teRXep4ldaM37frWtb9qxLfz8/8s/1rWt/tHH+q/WvJrvQ8HER8zWg/wCPr/tn/Wr61lQfaftX/LL/AFfv61fX7V/0x/Wvl8XuxUY6rUuLSXX/AB6v+H86hX7V/wBMf1pLn7V9lf8A1Pb19a+bxW572FWqNRamWqS/a/8Apj+tTL9r/wCmH618xikfQYYurUVp/wAeKfj/ADpi/a/+mH61Fafa/sSf6jv6+tfPYiOj1PocLuguOtZVxV64+15/5Y/rWVcfav8Apj+tTRj5n0WHZl3H/H4f+uf9aoXHWrVx9o+2H/Vf6v39aoT+f/0z/WvXpLY9nDsz7jpWPe/w/wC8K1J/Ox/yz/Wsi8835fufeHrXqUFqe1SloZt196s2XpV+68zdztrOl3/7NerSWheKlpsVJe9fU3ww/wCSbaJ/1wP/AKG1fLEu7npX1P8ADD/km2h/9cD/AOhtXqYXqfDZ89I+pk/EhN+q6KP+mFz/AOhQ1xtxCEtpWPACE5P0rvPHib9Y0cf9O9z/AOhQ1y15BiynPpG38q+ZzSpbMbf4f0N8uq8uB5fU9I1fxAmn2cdzaxx3aNFJOWWbavlxruYqQDk+g788jFJfazqNpqU1pHo/mq0LNZyfaVX7TKF3eXjHyd/mPHFatzZ2t4iLdW0M6owdRKgYK3qM9DUH9jaXkn+zbPLRCEnyF5jAAC9Pu4A46cCvsj5Ux7TxY9zf2tq2nNHuAFzIZciCQySRhOnznfEwyOOh6GulqommafE8Dx2Nsj24KwssKgxg9l446np61boAKKKKACiiigDyT9of/kQtP/7Csf8A6Kmr5uWvpH9of/kQtP8A+wrH/wCipq+blr67IP4L9f8AIwq7kq1J/CajWpP4TX2+G+E5mXbX7q1sW1Y9r91a2LalV2PNxJrW/atOy/j/AN81mW/atOy/j/3zXkV9zxK+zNu361rW/asm361rW/avJr7Hg4g0YP8Aj6/7Z/1q+tUIP+Pr/tn/AFq+tfLYvdio7omWkuv+PV/w/nSrSXX/AB6v+H86+bxW572E3RdWplqFamWvl8WfQ4YmWorT/jxT8f51KtRWn/Hin4/zr53EbM+hwu6K9x1rKuK1bjrWVcUqJ9HhzEuP+Pw/9c/61QuOtX7j/j8P/XP+tULjrXr0uh7WGMq46Vj3v8P+8K2LjpWPe/w/7wr1KG57NLYyrr71ZsvStK6+9WbL0r1aWxpitilL3r6m+GH/ACTbRP8Argf/AENq+WZe9fU3ww/5Jton/XA/+htXq4XqfC59tH1GeME363pQx0trn/0KGuf1CDGmXRx0hf8Aka6nxInma9pg/wCnW5/9ChrI1aHbo183PFvIen+ya+LzmpbNbf4f0OTC1eXD8vqd9RWLq+v/ANnWcdxDatMGiknZZS0JEca7m4K53dMAgfUVLearcWurxWK2iOs8TtFIZHGXVSdp+QjnB6EtjnbjNfoB4xq0Vz0HiWeWHT5nsFSO6uHtpP3j5iYOUHWMAZIHDFDzgAniuhoAKKKKACiqWryXUWlXD2YbzwvylF3MBnkqvOSBkgYOSOhrm7jUtf8AsNqbcXe5JZxIzWR8yba48lCNuEDoeXwApHODkUAcr+0P/wAiDp//AGFY/wD0VNXzapFfR/7Q5l/4QnThsTyf7TQlt53bvKl4xjpjPOfwr5xWvrcgv7F27/5GFXckUipNw2mmLUn8Jr7bD83LuczLlq67V5rWt5E9ayrX7q1sW1FXmtuediLGnbypx836VpWU8Y3/ADfxntVC37Vp2X8f++a8ive54le1ma1vcRZ+9+hrWt7mEY+f9DWdb9a1rftXk172PBxHKW4LqH7V9/8A5Z+h9avrdwf3/wBDVaD/AI+v+2f9avrXy+L3YqPLdCLeQf3/ANDSXN5B9lf5/TsfWrC0l1/x6v8Ah/Ovm8Vue9hbXRMt5B/z0/Q1Mt7b/wDPT9DSrUy18xirH0GGGre2/wDz0/8AHTUVpe2/2JP3nr2PrV1aitP+PFPx/nXz2I5bM+hwu6KFxeQE/wCs/Q1lXF3B/f8A0Nblx1rKuKmjy9j6LD3OfuLmH7Yfn/5Z+h9aoT3EX9/9DWpcf8fh/wCuf9aoXHWvXpW0PZw9zGnnjx979DWReTR/L838Q7Vt3HSse9/h/wB4V6lC1z2qXNYx7qVC3X9KzpXX1rUuvvVmy9K9WlaxeK5rFGV155r6n+GH/JNtD/64H/0Nq+WZe9fU3ww/5Jton/XA/wDobV6mF6nw2fXtH1LmtgHxBpuf+fW5/wDQ4az9ZC/2HqH/AF7Sf+gmr+unHiDTf+vW5/8AQ4azdYb/AIkl/wD9e0n/AKCa+MzmjzZrzf4f0PnlX5VynXXNna3iIt1bQzqjB1EqBgreoz0NQf2Rpgbd/Z1pnyvJz5C/6vGNvT7uABjpV2ivvySnHpOmxPC8en2iNBnySsKgx5JJ28cck9PWrlFFABRRRQAUUUUAeSftD/8AIhaf/wBhWP8A9FTV83LX0j+0P/yIWn/9hWP/ANFTV83LX12QfwX6/wCRhV3JVqT+E1GtSfwmvt8N8JzMu2v3VrYtqx7X7q1sW1KrsebiTWt+1adl/H/vmsy37Vp2X8f++a8ivueJX2Zt2/Wta37Vk2/Wta37V5NfY8HEGjB/x9f9s/61fWqEH/H1/wBs/wCtX1r5bF7sVHdEy0l1/wAer/h/OlWkuv8Aj1f8P5183itz3sJui6tTLUK1MtfL4s+hwxMtRWn/AB4p+P8AOpVqK0/48U/H+dfO4jZn0OF3RXuOtZVxWrcdayrilRPo8OYlx/x+H/rn/WqFx1q/cf8AH4f+uf8AWqFx1r16XQ9rDGVcdKx73+H/AHhWxcdKx73+H/eFepQ3PZpbGVdferNl6VpXX3qzZelerS2NMVsUpe9fU3ww/wCSbaJ/1wP/AKG1fLMvevqb4Yf8k20T/rgf/Q2r1cL1Phc+2j6ljxAca/pn/Xrc/wDoUNZervjRL8/9O8nT/dNaXiQ417TP+vW5/wDQoayNVb/iT3v/AF7yf+gmvDzChzY7m9D4jEVuWvy+hveINem02xhubURqJEkcfaomXeVXKxBSVIdjwM+h4NPfXEOttZx6jp4iksDdRMeSuCPmPz/MmDnjHA61u0V9Qemcxp3ic3aeG5HvNPI1S33SIh5MmzJ2Hd0DZGME+9dPRRQAUUUUAFFFFAHkn7Q//Ig6f/2FY/8A0VNXzaoFfSX7Q/8AyIWn/wDYVj/9FTV83LX1uQJOi79/8jCruSKBUm0bTTFqT+E19th4x5djmZctUXavFa1vGnpWVa/dWti2oqxjbY87ENmnbxJx8v61pWUMZ3/L/Ge9ULftWnZfx/75ryK8VfY8Su3Zmtb28Wfu/qa1re2hOPk/U1nW/Wta37V5NdKx4OIlLuW4LWH7V9z/AJZ+p9avraQf3P1NVoP+Pr/tn/Wr618vi92KjKV1qItnB/c/U0lzZwfZX+T07n1qwtJdf8er/h/Ovm8U9T3sK3dEy2cH/PP9TUy2Vv8A88/1NKtTLXzGKbPoMMxq2Vv/AM8//HjUVpZW/wBiT936/wAR9aurUVp/x4p+P86+exEpWep9Dhd0ULizgB/1f6msq4tIP7n6mty461lXFTRlLufRYdI5+4tofth+T/ln6n1qhPbxf3P1Nalx/wAfh/65/wBaoXHWvXpSemp7OHijGngjx939TWReQx/L8v8AEO9bdx0rHvf4f94V6lBu57VKMbbGPdRIG6frWdKi+lal196s2XpXq0m7F4qMbbFGVF54r6n+GH/JNtD/AOuB/wDQ2r5Zl719TfDD/km2if8AXA/+htXqYXqfDZ8klH1H+JzjXdL/AOva5/8AQoaxtUb/AIlN7/1wf/0E1reKzjXNL/69rn/0KGsPU2/4lN5wT+4fj/gJqK1Dmq8x+Z5hW5cdy+h6TRVOe9uIrWOZNLu5nc/NDG0QdPrucL+RNZOm+IZrnV5Le6EFumHBhfAlt2DoiLIQxUmTflcY6cbutd59SdFRXJt4qMemvcvqOl5h1VLORuisjSBcD5/lYKSc5I+U8emtpuqm71nVbB7i1ka0dNqxcOqsucMNx5z34+lAGtRRRQBU1K9/s/T5boR+YyYCpu2gsSAMnsMkZPasa88Xw2FjZy3EUST3F4bUxtOFVds3lO4YgZAJBAxk5HTkjoZYo54niljWSN1KujjIYHqCO4qGLTrKCAww2dvHEXVyiRKF3DGDgDqNq4PsPSgDyr9oeVR4J06HD7jqaNnYduPKlH3sYzz0zmvnFa+kf2h/+RB0/wD7Csf/AKKmr5tXPtX1uQP9y/X/ACMKu5MtSfwmolz7VJztPSvtsPL3djmZetfurWxbVjWu7avSta38z/Zoqy02POxCNm37Vp2X8f8Avmsm383j7n61pWXnfPjZ98+teRXep4ldaM37frWtb9qxLfz8/wDLP9a1rf7Rx/qv1rya70PBxEfM1oP+Pr/tn/Wr61lQfaftX/LL/V+/rV9ftX/TH9a+Xxe7FRjqtS4tJdf8er/h/OoV+1f9Mf1pLn7V9lf/AFPb19a+bxW572FWqNRamWqS/a/+mP61Mv2v/ph+tfMYpH0GGLq1Faf8eKfj/OmL9r/6YfrUVp9r+xJ/qO/r6189iI6PU+hwu6C461lXFXrj7Xn/AJY/rWVcfav+mP61NGPmfRYdmXcf8fh/65/1qhcdatXH2j7Yf9V/q/f1qhP5/wD0z/WvXpLY9nDsz7jpWPe/w/7wrUn87H/LP9ayLzzfl+594etepQWp7VKWhm3X3qzZelX7rzN3O2s6Xf8A7NerSWheKlpsVJe9fU3ww/5Jton/AFwP/obV8sS7uelfU/ww/wCSbaH/ANcD/wChtXqYXqfDZ89I+o3xcca3pX/Xtc/+hQ1g6i3/ABLLv/ri/wD6Ca2vGjhNa0kn/n3uf/Qoa5y/mB065HrE/wDI161PD88OY/Hs5quOa8v+H9D1aiiiuU++CiiigAooooAKKKKAPJP2h/8AkQtP/wCwrH/6Kmr5uWvpH9of/kQtP/7Csf8A6Kmr5uWvrsg/gv1/yMKu5KtSfwmo1qT+E19vhvhOZl21+6tbFtWPa/dWti2pVdjzcSa1v2rTsv4/981mW/atOy/j/wB815Ffc8Svszbt+ta1v2rJt+ta1v2rya+x4OINGD/j6/7Z/wBavrVCD/j6/wC2f9avrXy2L3YqO6JlpLr/AI9X/D+dKtJdf8er/h/Ovm8Vue9hN0XVqZahWplr5fFn0OGJlqK0/wCPFPx/nUq1Faf8eKfj/OvncRsz6HC7or3HWsq4rVuOtZVxSon0eHMS4/4/D/1z/rVC461fuP8Aj8P/AFz/AK1QuOtevS6HtYYyrjpWPe/w/wC8K2LjpWPe/wAP+8K9ShuezS2Mq6+9WbL0rSuvvVmy9K9WlsaYrYpS96+pvhh/yTbRP+uB/wDQ2r5Zl719TfDD/km2if8AXA/+htXq4XqfC59tH1M/4gyeXq2jn1guf/Qoa5G6uM2c4AJ/dtx+FdJ8TpPL1PRDnrDdf+hQ1xMlxuicZ6qRX2GAp3wd/U/IM4pc2ac3+H9D3GefUDaxyWdlEZmPzxXVx5ewfVFcE/5zWTpt9qx1eRdQSdY8OJIkhLRRNvRYvLcIC4KlixOcY524roqK+bPuzjLnUNcXTtTa3l1FpUn/ANDP2T5nG37pHkDjPqAOg8zrjZsLi+fxHfxO1w9iqAxmWLaqvk5AOxcj6F+mcjodqigAooooAKKiubmGztpLid9kUYyxwT+QHJPsOtU21ywS1iuTJL5MrlA4t5CFYNtIf5fkw3B3Y6H0oA81/aH/AORB0/8A7Csf/oqavm1T7GvpL9of/kQtP/7Csf8A6Kmr5uWvrcgv7F+v+RhV3JFPsakz8p4NMWpP4TX22HT5dzmZctW+VflP5VrW7/7LflWVa/dWti2oqp23POxBp28nT5H/ACrSspfv/u3++e1ULftWnZfx/wC+a8iunfc8Su1Zmtbzc/6uT/vmta3n6fupf++azrfrWtb9q8munY8HENdi3Bcf6V/qZf8AV/3fer63P/TGb/vmq0H/AB9f9s/61fWvl8XuxUWrrQRbn/pjN/3xSXNz/or/ALmbt/B71YWkuv8Aj1f8P5183itz3sK1dEy3P/TCb/vipluv+mE//fFKtTLXzGKaPoMMNW6/6YT/APfFRWl1/oSfuJ+/8HvV1aitP+PFPx/nXz2Ias9D6HC7ooXFzz/qJv8Avisq4uP+mM3/AHzW5cdayripotdj6LDnP3E/+mH91L/q/wC771Qnm/6ZSf8AfNalx/x+H/rn/WqFx1r16TWmh7OHTMaebj/Vyf8AfNZF5L935H+8O1bdx0rHvf4f94V6lBq57VJO25j3Unzfcb8qzpX/ANlvyrUuvvVmy9K9WlsXik7blGVuvB/Kvqf4Yf8AJNtD/wCuB/8AQ2r5Zl719TfDD/km2if9cD/6G1ephep8Nn20fU5z4tyeXqGhH1iuv5w158LjcQPXiu3+NEnl3vh8+sd1/OGvM4rgmZAOSWGM19/lVO+XX/xfqfmeYUubHc3ofU1FU55dTW1ja3tLSS4J/eRyXTIij2YRkn8hWTpt9qx1eRdQSdY8OJIkhLRRNvRYvLcIC4KlixOcY524r4w+pOioriDdeIUsHK6jfTyw3GY2OlFDdrsH7raU/djd/wAtDxz7GtjQrnV5tU1JNREnkiRvJBg8tY1EjqoVsfPlArZycZxxQBv0UUUAV76zjv7OS2lZlVwPmQ/MpByCM9wQD+FZn/COJ5UEX9p33lw3LXWw+ViR2bcdw2cjcSwHYnjouNuigDyH9oeJT4J06bL7hqaLjeduPKlP3c4zx1xmvnFa+kf2h/8AkQtP/wCwrH/6Kmr5uWvrsg/gv1/yMKu5KtSfwmo1qT+E19vhvhOZl21+6tbFtWPa/dWti2pVdjzcSa1v2rTsv4/981mW/atOy/j/AN815Ffc8Svszbt+ta1v2rJt+ta1v2rya+x4OINGD/j6/wC2f9avrVCD/j6/7Z/1q+tfLYvdio7omWkuv+PV/wAP50q0l1/x6v8Ah/Ovm8Vue9hN0XVqZahWplr5fFn0OGJlqK0/48U/H+dSrUVp/wAeKfj/ADr53EbM+hwu6K9x1rKuK1bjrWVcUqJ9HhzEuP8Aj8P/AFz/AK1QuOtX7j/j8P8A1z/rVC46169Loe1hjKuOlY97/D/vCti46Vj3v8P+8K9ShuezS2Mq6+9WbL0rSuvvVmy9K9WlsaYrYpS96+pvhh/yTbRP+uB/9DavlmXvX1N8MP8Akm2if9cD/wChtXq4XqfC59tH1OI+Osnl3Ph056pdfzhryq0uN17AM9ZFH616X+0FJ5c3hs+q3X/tGvINOnzqdoPWZP5iv07Jad8pv5S/U+FxNLmr83ofZtFFFfnp6wUUUUAFFFFABRRRQB5J+0P/AMiDp/8A2FY//RU1fNqgV9JftD/8iFp//YVj/wDRU1fNy19bkCTou/f/ACMKu5IoFSbRtNMWpP4TX22HjHl2OZly1Rdq8VrW8aelZVr91a2LairGNtjzsQ2advEnHy/rWlZQxnf8v8Z71Qt+1adl/H/vmvIrxV9jxK7dma1vbxZ+7+prWt7aE4+T9TWdb9a1rftXk10rHg4iUu5bgtYftX3P+WfqfWr62kH9z9TVaD/j6/7Z/wBavrXy+L3YqMpXWoi2cH9z9TSXNnB9lf5PTufWrC0l1/x6v+H86+bxT1Pewrd0TLZwf88/1NTLZW//ADz/AFNKtTLXzGKbPoMMxq2Vv/zz/wDHjUVpZW/2JP3fr/EfWrq1Faf8eKfj/OvnsRKVnqfQ4XdFC4s4Af8AV/qayri0g/ufqa3LjrWVcVNGUu59Fh0jn7i2h+2H5P8Aln6n1qhPbxf3P1Nalx/x+H/rn/WqFx1r16Unpqezh4oxp4I8fd/U1kXkMfy/L/EO9bdx0rHvf4f94V6lBu57VKMbbGPdRIG6frWdKi+lal196s2XpXq0m7F4qMbbFGVF54r6n+GH/JNtD/64H/0Nq+WZe9fU3ww/5Jton/XA/wDobV6mF6nw2fJJR9Tzn9oo4k8M/S7/APaNeM6Yx/tWzxyfPTGT/tCvZP2jDiTwz9Lv/wBo14xpTf8AE3sv+u8f/oQr9OyWty5Ty/4v1Pkai/eH2hOupy2sf2eW0tbjOZPMiadMegwyH8f0rJ02DWbbV5JbwTzKQ6yukgEUhZ08to4y52BU37hwT/tnFdFRX56dZy0uh3stheRrLqEb3N6BHjUJd0EIIQsG391DOB6uARxxPo1rq8PiDUZb3z/skhfyt825fvnbgbz/AA+ipjod5+auiooAKKKKACiqWrx3UulXEdmWE5X5QjbWIzyFbjBIyAcjBPUVzs1pr0mnRRxx3yOhuVjzdruVmYGB3YP86IpKsCWJI6NwaAOT/aH/AORC0/8A7Csf/oqavm5a+u/iF4Kk8dadY6Y96La0iuhPKyx7nyI3UYJOMZYDGPxrhh+zppw/5mG7/wDAdf8AGvfyrMaGFp8tS97mU4OT0PAlqT+E170P2d9OH/MwXX/gOv8AjTv+GedPx/yMF1/4Dr/jX01HibAQVm39xi6MjxC1+6tbFtXrkfwBsowANfuOPW3X/GpJPgfHDbyNb63JJKEJRGgADNjgE7uOaKnEuAktG/uOOtgqs9jzW37Vp2X8f++a9ET4NRp012T8bUf/ABVTRfCQQ52643Jzzaj/AOKrz6ue4ST0b+48yrlGJktLfechb9a1rftXRp8M5E6a0PxtP/s6sJ8P7lOmsp+Np/8AZ1wVc1w8tr/ceXW4dxstkvvMWD/j6/7Z/wBavrWingi8STeNYizjH/Hmf/jlOi8J6oXlD6pbqqviMi0J3LtHJ/ecc5H4V4eIqxnflFT4dxsXql95SWkuv+PV/wAP51qDwpqA/wCYvB/4BH/45RJ4T1CSModXgwfSyP8A8crxq+HnPY9WhlGJg1e33kC1MtTDw3qQ/wCYtbf+AR/+OU4eHtTH/MWtv/AI/wDxyvEr5RiZ7W+89ajhKkNxq1Faf8eKfj/OpToesi4RV1K0MRRiz/YzkNkYGPM7gt+XvTo/D2qRRCNdWtsD1sj/APHK8etw3jp7JfeevQkoNXKFx1rKuK6F/DWpP11a2/8AAI//AByoH8HXz9dYh/CzP/xyinw3jo7pfeexSx9GO5xVx/x+H/rn/WqFx1run8AXLybzrMecY/48/wD7Oqtz8ObsoDHq6O29QQbXHylhuP3+wyfwrvhkeLW6X3no0c4w0N7/AHHnVx0rHvf4f94V60/wsZ+ut/8Akr/9nVaX4PpLjdrj8HPFqP8A4qu6lleIi9bfeelT4iwUVq39x4zdferNl6V7hJ8E4JDltel/C2H/AMVULfAq1brr8/8A4DD/AOKruhg6sVqVX4jwM1o39x4RL3r6m+GH/JNtE/64H/0Nq4tvgFZN11+4/wDAdf8AGvR/CejP4e8M2mkySCT7LvRXxjcu9ipI9dpGfeu+hSlC9z5jNcbSxKXs+h5F+0ecP4Y+l3/7RrxbSW/4nNj/ANfEf/oQr6q+IPw4tviAdONxqE1n9h83b5cYbdv2Zzn02D86463/AGebC2uYp4/EFyXicOu63UjIORnmvp8HmdOjhfYy31/E8GUG5XPZ6Kpzxam1rGtvd2kdwD+8kktWdGHsokBH5mqFguo2+s3xuLe5lheOALIJF8tpMkOyI0hKDBU4/wBk4yevgGpt0Vwzab4lOn7B9u8zywGAvQGa48tgZQ2/iLeVOzjp93tW5pNtqkWr3Ml2Z/KPm7mebckmZMxeWuTs2pweFycdetAG7RRRQAUUUUAcjD4av4or2HbZFJxFu/eN/pDJKXJkGz/lorFW64wPvA8dBo9lJp+mR20uwMrOwSM5SMM5YIvA4UEKOBwOg6VeooAKKKKACsC50a7bxPDqcS2zosisXdysgTYyNGMKfly2/rywxjuN+igDG8P6ZPpy3j3FtZ2z3E3mCOzcmMDGBxsXnAGTzk+gwBs0UUAFFFFAGbrenyajZRxJHBMEmWRoLg4jlA/hY4P16HkCsiPw9qhn8Oyy3lsRpYCvGys+7Ebxl1bI+Zgy8EfLzgnnPU0UAFFFFABUV1E09pNCkhieRGVZF6qSMZH0qWigDjf+EVvX0p7GS3sltzcpP9kiuHSNwIthQkJlRuAkyAckke56uxhlt7C2hnkEs0cSo8gXaGYAAkDtk9qnooAKKKKACuXfQL8apfXEQtNtzFcRmZpGDv5gUpuAHOwptHzdGJBBGD1FFAGL4Y0ifRdMe2n8oEyl1ERBABAH8KIucg9FUfU5J2qKKACiiigDJ1/T7jULW2W1gtZZobuGcNcOU2BJFYlSFY5IBHbr1qtFo14fEKanKtspLiR3Ry0gHk7DAPlGY93z5yOf4e9b9FABRRRQAUUUUAFFFFABRRRQB//Z", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "op.setConstraints(rhs)\n", "op.setConstraints(sol)\n", "rk = sol.copy(\"residual\")\n", "def cb(xk): # a callback to print the residual norm in each step\n", " x_h = space.function(\"iterate\", dofVector=xk)\n", " op(x_h,rk)\n", " print(rk.as_numpy[:].dot(rk.as_numpy[:]), flush=True, end='\\t')\n", "sol.as_numpy[:], _ = solver(lin.as_numpy, rhs.as_numpy, x0=sol.as_numpy,\n", " callback=cb, tol=1e-10, atol=1e-10)\n", "sol.plot()" ] }, { "cell_type": "markdown", "id": "11a15922", "metadata": {}, "source": [ ".. index::\n", " pair: Solvers; Petsc\n", "\n", "## Using PETSc and Petsc4Py\n", "The following requires that a PETSc installation was found during the\n", "configuration of ``dune``. Furthermore some examples make use of the\n", "Python package ``petsc4py`." ] }, { "cell_type": "code", "execution_count": 13, "id": "d631452d", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:13:21.599909Z", "iopub.status.busy": "2024-03-01T11:13:21.599665Z", "iopub.status.idle": "2024-03-01T11:13:21.621010Z", "shell.execute_reply": "2024-03-01T11:13:21.619607Z" } }, "outputs": [], "source": [ "from dune.common.checkconfiguration import assertCMakeHave, ConfigurationError\n", "try:\n", " assertCMakeHave(\"HAVE_PETSC\")\n", " petsc = True\n", "except ConfigurationError:\n", " print(\"Dune not configured with petsc - skipping example\")\n", " petsc = False\n", "try:\n", " import petsc4py\n", " petsc4py.init(sys.argv)\n", " from petsc4py import PETSc\n", "except ModuleNotFoundError:\n", " print(\"petsc4py module not found -- skipping example\")\n", " petsc4py = None" ] }, { "cell_type": "markdown", "id": "d62eb52f", "metadata": {}, "source": [ "Switching to a storage based on the PETSc solver package and solving the\n", "system using the dune-fem bindings can be achieved by using the\n", "``storage`` argument to the space constructor" ] }, { "cell_type": "code", "execution_count": 14, "id": "0b9529f1", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:13:21.626481Z", "iopub.status.busy": "2024-03-01T11:13:21.626145Z", "iopub.status.idle": "2024-03-01T11:15:18.070098Z", "shell.execute_reply": "2024-03-01T11:15:18.068672Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer(petsc) size: 64 L^2, H^1 error: 1.93079e-05, 9.99164e-04\n" ] } ], "source": [ "if petsc:\n", " spacePetsc = solutionSpace(gridView, order=2, storage='petsc')\n", " # first we will use the petsc solver available in the `dune-fem` package\n", " # (using the sor preconditioner)\n", " schemePetsc = solutionScheme(a == b, space=spacePetsc,\n", " parameters={\"linear.preconditioning.method\":\"sor\"})\n", " dt.value = scheme.model.dt\n", " u_h = spacePetsc.interpolate(initial, name='u_h')\n", " u_h_n = u_h.copy(name=\"previous\")\n", " evolve(schemePetsc, u_h, u_h_n, endTime)\n", " error = u_h - exact_end\n", " print(\"Forchheimer(petsc) size: \", gridView.size(0), \"L^2, H^1 error:\",'{:0.5e}, {:0.5e}'.format(\n", " *[ sqrt(e) for e in integrate([error**2,inner(grad(error),grad(error))]) ]))" ] }, { "cell_type": "markdown", "id": "5e69a030", "metadata": { "lines_to_next_cell": 0 }, "source": [ "Implementing a Newton Krylov solver using the binding provided by petsc4py" ] }, { "cell_type": "code", "execution_count": 15, "id": "76082191", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:15:18.075761Z", "iopub.status.busy": "2024-03-01T11:15:18.075536Z", "iopub.status.idle": "2024-03-01T11:15:37.918316Z", "shell.execute_reply": "2024-03-01T11:15:37.916872Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer(petsc) size: 64 L^2, H^1 error: 1.93091e-05, 9.99162e-04\n" ] } ], "source": [ "if petsc4py is not None and petsc:\n", " class Scheme3:\n", " def __init__(self, scheme):\n", " self.model = scheme.model\n", " self.jacobian = scheme.linear()\n", " self.ksp = PETSc.KSP()\n", " self.ksp.create(PETSc.COMM_WORLD)\n", " # use conjugate gradients method\n", " self.ksp.setType(\"cg\")\n", " # and incomplete Cholesky\n", " self.ksp.getPC().setType(\"icc\")\n", " self.ksp.setOperators(self.jacobian.as_petsc)\n", " self.ksp.setFromOptions()\n", " def solve(self, target):\n", " res = target.copy(name=\"residual\")\n", " sol_coeff = target.as_petsc\n", " res_coeff = res.as_petsc\n", " n = 0\n", " while True:\n", " schemePetsc(target, res)\n", " absF = numpy.sqrt( res_coeff.dot(res_coeff) )\n", " if absF < 1e-10:\n", " break\n", " schemePetsc.jacobian(target, self.jacobian)\n", " self.ksp.solve(res_coeff, res_coeff)\n", " sol_coeff -= res_coeff\n", " n += 1\n", "\n", " u_h.interpolate(initial)\n", " scheme3_cls = Scheme3(schemePetsc)\n", " evolve(scheme3_cls, u_h, u_h_n, endTime)\n", " error = u_h - exact_end\n", " print(\"Forchheimer(petsc) size: \", gridView.size(0), \"L^2, H^1 error:\",'{:0.5e}, {:0.5e}'.format(\n", " *[ sqrt(e) for e in integrate([error**2,inner(grad(error),grad(error))]) ]))" ] }, { "cell_type": "markdown", "id": "207908a4", "metadata": { "lines_to_next_cell": 0 }, "source": [ "Using the petsc4py bindings for the non linear KSP solvers from PETSc" ] }, { "cell_type": "code", "execution_count": 16, "id": "ae517044", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:15:37.924116Z", "iopub.status.busy": "2024-03-01T11:15:37.923858Z", "iopub.status.idle": "2024-03-01T11:16:13.637055Z", "shell.execute_reply": "2024-03-01T11:16:13.635699Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Forchheimer(petsc4py) size: 64 L^2, H^1 error: 1.93091e-05, 9.99162e-04\n" ] } ], "source": [ "if petsc4py is not None and petsc is not None:\n", " class Scheme4:\n", " def __init__(self, scheme):\n", " self.model = scheme.model\n", " self.res = scheme.space.interpolate([0],name=\"residual\")\n", " self.scheme = scheme\n", " self.jacobian = scheme.linear()\n", " self.snes = PETSc.SNES().create()\n", " self.snes.setFunction(self.f, self.res.as_petsc.duplicate())\n", " self.snes.setUseMF(False)\n", " self.snes.setJacobian(self.Df, self.jacobian.as_petsc, self.jacobian.as_petsc)\n", " self.snes.getKSP().setType(\"cg\")\n", " self.snes.setFromOptions()\n", "\n", " def f(self, snes, x, f):\n", " # setup discrete function using the provide petsc vectors\n", " inDF = self.scheme.space.function(\"tmp\",dofVector=x)\n", " outDF = self.scheme.space.function(\"tmp\",dofVector=f)\n", " self.scheme(inDF,outDF)\n", "\n", " def Df(self, snes, x, m, b):\n", " inDF = self.scheme.space.function(\"tmp\",dofVector=x)\n", " self.scheme.jacobian(inDF, self.jacobian)\n", " return PETSc.Mat.Structure.SAME_NONZERO_PATTERN\n", "\n", " def solve(self, target):\n", " sol_coeff = target.as_petsc\n", " self.res.clear()\n", " self.snes.solve(self.res.as_petsc, sol_coeff)\n", "\n", " u_h.interpolate(initial)\n", " scheme4_cls = Scheme4(schemePetsc)\n", " evolve(scheme4_cls, u_h, u_h_n, endTime)\n", " error = u_h - exact_end\n", " print(\"Forchheimer(petsc4py) size: \", gridView.size(0), \"L^2, H^1 error:\",'{:0.5e}, {:0.5e}'.format(\n", " *[ sqrt(e) for e in integrate([error**2,inner(grad(error),grad(error))]) ]))" ] }, { "cell_type": "markdown", "id": "fd41e660", "metadata": {}, "source": [ ".. index::\n", " triple: I/O; Logging; Parameters\n", "\n", "## Accessing and reusing values of parameters\n", "Sometimes it is necessary to extract which parameters were read and which\n", "values were used, e.g., for debugging purposes like finding spelling\n", "in the parameters provided to a scheme.\n", "Note that this information can only be reliably obtained after usage of\n", "the scheme, e.g., after calling solve as shown in the example below.\n", "To add logging to a set of parameters passed to a `scheme` one simply\n", "needs to add a `logging` key to the parameter dictionary provided to the scheme\n", "with a tag (string) that is used in the output.\n", "\n", "As an example we will solve the simple Laplace equation from the\n", "introduction but pass some preconditioning parameters to the scheme." ] }, { "cell_type": "code", "execution_count": 17, "id": "52392cdb", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:16:13.642369Z", "iopub.status.busy": "2024-03-01T11:16:13.642159Z", "iopub.status.idle": "2024-03-01T11:19:00.869390Z", "shell.execute_reply": "2024-03-01T11:19:00.867965Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== Dune::CGSolver\n", " Iter Defect Rate\n", " 0 0.19886\n", " 1 0.138041 0.694164\n", " 2 0.0533951 0.386805\n", " 3 0.0335279 0.627921\n", " 4 0.00936836 0.27942\n", " 5 0.0158677 1.69376\n", " 6 0.00553746 0.348977\n", " 7 0.00407438 0.735784\n", " 8 0.00238342 0.584978\n", " 9 0.00097925 0.410859\n", " 10 0.00030568 0.312157\n", " 11 0.000177396 0.580333\n", " 12 4.54516e-05 0.256215\n", " 13 2.72129e-05 0.598721\n", " 14 2.549e-05 0.936688\n", " 15 1.32301e-05 0.51903\n", " 16 7.51142e-06 0.567755\n", " 17 2.09346e-06 0.278703\n", " 18 1.48491e-06 0.70931\n", " 19 7.68912e-07 0.517818\n", " 20 2.79925e-07 0.364054\n", " 21 5.71861e-08 0.20429\n", " 22 2.5419e-08 0.444496\n", " 23 1.38267e-08 0.543951\n", " 24 6.60951e-09 0.478026\n", " 25 1.93829e-09 0.293257\n", " 26 6.68847e-10 0.345072\n", " 27 3.43101e-10 0.512974\n", " 28 1.06162e-10 0.309418\n", " 29 2.61402e-11 0.24623\n", " 30 5.42124e-12 0.207391\n", " 31 5.98748e-12 1.10445\n", " 32 2.10778e-12 0.352031\n", " 33 9.47556e-13 0.449552\n", "=== rate=0.453848, T=0.237529, TIT=0.00719784, IT=33\n" ] } ], "source": [ "import dune.fem\n", "from dune.grid import structuredGrid\n", "from dune.fem.space import lagrange\n", "from dune.fem.scheme import galerkin\n", "from ufl import (TestFunction, TrialFunction, SpatialCoordinate,\n", " dx, grad, inner, dot, sin, cos, pi )\n", "gridView = structuredGrid([0, 0], [1, 1], [200, 200])\n", "space = lagrange(gridView, order=1, storage=\"istl\")\n", "u_h = space.interpolate(0, name='u_h')\n", "x = SpatialCoordinate(space)\n", "u = TrialFunction(space)\n", "v = TestFunction(space)\n", "\n", "f = (8*pi**2+1) * cos(2*pi*x[0])*cos(2*pi*x[1])\n", "a = ( inner(grad(u),grad(v)) + u*v ) * dx\n", "l = f*v * dx\n", "scheme = galerkin( a==l, solver=\"cg\", parameters=\n", " {\"newton.linear.tolerance\": 1e-12,\n", " \"newton.linear.verbose\": True,\n", " \"newton.linear.preconditioning.method\": \"amg-ilu\",\n", " \"fem.solver.newton.linear.errormeasure\": \"relative\",\n", " \"logging\": \"precon-amg\"\n", " } )\n", "info = scheme.solve(target=u_h)" ] }, { "cell_type": "markdown", "id": "3faf98f3", "metadata": {}, "source": [ "We use the `pprint` (pretty print) module if available to get nicer\n", "output." ] }, { "cell_type": "code", "execution_count": 18, "id": "952a26aa", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:19:00.874841Z", "iopub.status.busy": "2024-03-01T11:19:00.874618Z", "iopub.status.idle": "2024-03-01T11:19:00.881517Z", "shell.execute_reply": "2024-03-01T11:19:00.880261Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'default': {('fem.threads.communicationthread', 'false'), ('fem.dofmanager.clearresizedarrays', 'true'), ('fem.dofmanager.memoryfactor', '1.1')},\n", " 'precon-amg': {('fem.solver.newton.lineSearch', 'none'),\n", " ('fem.solver.newton.linear.errormeasure', 'absolute'),\n", " ('fem.solver.newton.linear.matrix.overflowfraction', '1'),\n", " ('fem.solver.newton.linear.maxiterations', '2147483647'),\n", " ('fem.solver.newton.linear.preconditioning.iterations', '1'),\n", " ('fem.solver.newton.linear.preconditioning.relaxation', '1.1'),\n", " ('fem.solver.newton.linear.threading', 'true'),\n", " ('fem.solver.newton.linear.tolerance.strategy', 'none'),\n", " ('fem.solver.newton.maxiterations', '2147483647'),\n", " ('fem.solver.newton.maxlinesearchiterations', '2147483647'),\n", " ('fem.solver.newton.tolerance', '1e-06'),\n", " ('fem.solver.newton.verbose', 'false'),\n", " ('newton.linear.method', 'cg'),\n", " ('newton.linear.preconditioning.method', 'amg-ilu'),\n", " ('newton.linear.tolerance', '1e-12'),\n", " ('newton.linear.verbose', 'True')},\n", " 'program code': {('fem.adaptation.method', 'callback')}}\n" ] } ], "source": [ "try:\n", " from pprint import pprint as _pprint\n", " pprint = lambda *args,**kwargs: _pprint(*args,**kwargs,width=200,compact=False)\n", "except ImportError:\n", " pprint = print\n", "\n", "pprint(dune.fem.parameter.log())" ] }, { "cell_type": "markdown", "id": "1a91f78f", "metadata": {}, "source": [ "Note above that all parameters are printed including some default ones\n", "used in other parts of the code. If multiple schemes with different\n", "`logging` parameter strings are used, all would be shown using the `log`\n", "method as shown above.\n", "To access only the parameters used in the scheme simply use\n", "either `dune.fem.parameter.log()[\"tag\"])` or access the parameter log\n", "through the scheme:" ] }, { "cell_type": "code", "execution_count": 19, "id": "ac03c675", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:19:00.886572Z", "iopub.status.busy": "2024-03-01T11:19:00.886330Z", "iopub.status.idle": "2024-03-01T11:19:00.892311Z", "shell.execute_reply": "2024-03-01T11:19:00.891240Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{('fem.solver.newton.lineSearch', 'none'),\n", " ('fem.solver.newton.linear.errormeasure', 'absolute'),\n", " ('fem.solver.newton.linear.matrix.overflowfraction', '1'),\n", " ('fem.solver.newton.linear.maxiterations', '2147483647'),\n", " ('fem.solver.newton.linear.preconditioning.iterations', '1'),\n", " ('fem.solver.newton.linear.preconditioning.relaxation', '1.1'),\n", " ('fem.solver.newton.linear.threading', 'true'),\n", " ('fem.solver.newton.linear.tolerance.strategy', 'none'),\n", " ('fem.solver.newton.maxiterations', '2147483647'),\n", " ('fem.solver.newton.maxlinesearchiterations', '2147483647'),\n", " ('fem.solver.newton.tolerance', '1e-06'),\n", " ('fem.solver.newton.verbose', 'false'),\n", " ('newton.linear.method', 'cg'),\n", " ('newton.linear.preconditioning.method', 'amg-ilu'),\n", " ('newton.linear.tolerance', '1e-12'),\n", " ('newton.linear.verbose', 'True')}\n" ] } ], "source": [ "pprint(scheme.parameterLog())" ] }, { "cell_type": "markdown", "id": "0591427b", "metadata": {}, "source": [ "One can easily reuse these parameters to construct another scheme by\n", "converting the result of the above call to a dictionary.\n", "As an example change the above problem to a PDE with Dirichlet conditions\n", "but turn of verbose output of the solver.\n", "\n", ".. note:: the `logging` parameter has to be set if we want to use the\n", "`parameterLog` method on the scheme." ] }, { "cell_type": "code", "execution_count": 20, "id": "86b87556", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:19:00.897494Z", "iopub.status.busy": "2024-03-01T11:19:00.897286Z", "iopub.status.idle": "2024-03-01T11:20:06.633770Z", "shell.execute_reply": "2024-03-01T11:20:06.632168Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{('fem.solver.newton.lineSearch', 'none'),\n", " ('fem.solver.newton.linear.errormeasure', 'absolute'),\n", " ('fem.solver.newton.linear.matrix.overflowfraction', '1'),\n", " ('fem.solver.newton.linear.maxiterations', '2147483647'),\n", " ('fem.solver.newton.linear.preconditioning.iterations', '1'),\n", " ('fem.solver.newton.linear.preconditioning.relaxation', '1.1'),\n", " ('fem.solver.newton.linear.threading', 'true'),\n", " ('fem.solver.newton.linear.tolerance.strategy', 'none'),\n", " ('fem.solver.newton.maxiterations', '2147483647'),\n", " ('fem.solver.newton.maxlinesearchiterations', '2147483647'),\n", " ('fem.solver.newton.tolerance', '1e-06'),\n", " ('fem.solver.newton.verbose', 'false'),\n", " ('newton.linear.method', 'cg'),\n", " ('newton.linear.preconditioning.method', 'amg-ilu'),\n", " ('newton.linear.tolerance', '1e-12'),\n", " ('newton.linear.verbose', 'False')}\n" ] } ], "source": [ "param = dict(scheme.parameterLog()) # this method returns a set of pairs which we can convert to a dictionary\n", "param[\"logging\"] = \"Dirichlet\" # only needed to use the `parameterLog` method\n", "param[\"newton.linear.verbose\"] = False\n", "scheme2 = galerkin( [a==l,DirichletBC(space,0)], parameters=param )\n", "u_h.clear()\n", "info = scheme2.solve(target=u_h)\n", "pprint(scheme2.parameterLog())" ] }, { "cell_type": "markdown", "id": "14d97dfc", "metadata": {}, "source": [ "### Parameter hints\n", "\n", ".. tip:: To get information about available values for some parameters\n", "(those with string arguments) a possible approach is to provide a non valid\n", "string, e.g., `\"help\"`." ] }, { "cell_type": "code", "execution_count": 21, "id": "4916f197", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:20:06.639293Z", "iopub.status.busy": "2024-03-01T11:20:06.639066Z", "iopub.status.idle": "2024-03-01T11:20:06.791992Z", "shell.execute_reply": "2024-03-01T11:20:06.790530Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help for parameter 'fem.solver.newton.linear.preconditioning.method':\n", "Valid values are: none, ssor, sor, ilu, gauss-seidel, jacobi, amg-ilu, amg-jacobi, ildl\n", "\n" ] } ], "source": [ "scheme = galerkin( a==l, solver=\"cg\", parameters=\n", " {\"newton.linear.tolerance\": 1e-12,\n", " \"newton.linear.verbose\": True,\n", " \"newton.linear.preconditioning.method\": \"help\",\n", " \"fem.solver.newton.linear.errormeasure\": \"relative\",\n", " \"logging\": \"precon-amg\"\n", " } )\n", "try:\n", " scheme.solve(target=u_h)\n", "except RuntimeError as rte:\n", " print(rte)" ] }, { "cell_type": "markdown", "id": "bd0bac33", "metadata": {}, "source": [ ".. index::\n", " triple: Solvers; Available solvers; Parameters\n", "\n", "## Available solvers and parameters\n", "Upon creation of a discrete function space one also have to specifies the\n", "storage which is tied to the solver backend.\n", "As mentioned, different linear algebra backends can be accessed by changing setting the\n", "`storage` parameter during construction of the discrete space. All\n", "discrete functions and operators/schemes based on this space will then\n", "use this backend. Available backends are `numpy,istl,petsc`.\n", "Note that not all methods which are available in `dune-istl` or `PETSc` have been forwarded\n", "to be used with `dune-fem`." ] }, { "cell_type": "code", "execution_count": 22, "id": "8a6ce5a1", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:20:06.797062Z", "iopub.status.busy": "2024-03-01T11:20:06.796852Z", "iopub.status.idle": "2024-03-01T11:20:06.802544Z", "shell.execute_reply": "2024-03-01T11:20:06.801314Z" } }, "outputs": [], "source": [ "space = solutionSpace(gridView, order=2, storage='numpy')" ] }, { "cell_type": "markdown", "id": "2dd9ca63", "metadata": {}, "source": [ "Switching is as simple as passing `storage='istl'` or `storage='petsc'`.\n", "Here is a summary of the available backends\n", "\n", "| Solver | Description |\n", "| --------- | ------------------------------------------------------------------- |\n", "| numpy | the storage is based on a raw C pointer which can be |\n", "| | directly accessed as a numpy.array using the Python buffer protocol |\n", "| | To change the underlying vector of a discrete function 'u_h' use |\n", "| | 'uh.as_numpy[:]'. |\n", "| | As shown in the examples, linear operators return a |\n", "| | `scipy.sparse.csr_matrix` through the 'as_numpy' property. |\n", "| istl | data is stored in a block vector/matrix from the dune.istl package |\n", "| | Access through 'as_istl' |\n", "| petsc | data is stored in a petsc vector/matrix which can also be used with |\n", "| | petsc4py on the python side using 'as_petsc' |" ] }, { "cell_type": "markdown", "id": "a3459bf4", "metadata": {}, "source": [ ".. index:: Solvers; Parameters\n", "\n", "When creating a scheme, there is the possibility to select a linear\n", "solver for the internal Newton method.\n", "In addition the behavior of the solver can be customized through a\n", "parameter dictionary. This allows to set tolerances, verbosity, but also\n", "which preconditioner to use.\n", "\n", "For details see the help available for a scheme:" ] }, { "cell_type": "code", "execution_count": 23, "id": "5c2f2dad", "metadata": { "execution": { "iopub.execute_input": "2024-03-01T11:20:06.807506Z", "iopub.status.busy": "2024-03-01T11:20:06.807307Z", "iopub.status.idle": "2024-03-01T11:20:06.817467Z", "shell.execute_reply": "2024-03-01T11:20:06.816180Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on Scheme in module dune.generated.femscheme_9b7005253a7d0b23b665ce4b07198dde object:\n", "\n", "class Scheme(pybind11_builtins.pybind11_object)\n", " | A scheme finds a solution `u=ufl.TrialFunction` for a given variational equation.\n", " | The main method is `solve` which takes a discrete functions as `target` argument to\n", " | store the solution. The method always uses a Newton method to solve the problem.\n", " | The linear solver used in each iteration of the Newton method can be chosen\n", " | using the `solver` parameter in the constructor of the scheme. Available solvers are:\n", " | ------------------------------------------\n", " | | Solver | Storage |\n", " | | name | numpy | istl | petsc |\n", " | |----------|---------|---------|---------|\n", " | | bicg | --- | --- | x |\n", " | | bicgstab | x | x | x |\n", " | | cg | x | x | x |\n", " | | gmres | x | x | x |\n", " | | gradient | --- | x | --- |\n", " | | loop | --- | x | --- |\n", " | | minres | --- | x | x |\n", " | | preonly | --- | --- | x |\n", " | | superlu | --- | x | --- |\n", " | ------------------------------------------\n", " | \n", " | In addition the direct solvers from the `suitesparse` package can be used with the\n", " | `numpy` storage. In this case provide a tuple as `solver` argument with \"suitesparse\" as\n", " | first argument and the solver to use as second, e.g.,\n", " | 'solver=(\"suitesparse\",\"umfpack\")'.\n", " | \n", " | The detailed behavior of the schemes can be customized by providing a\n", " | `parameters` dictionary to the scheme constructor, e.g.,\n", " | {\"newton.tolerance\": 1e-3, # tolerance for newton solver\n", " | \"newton.verbose\": False, # toggle iteration output\n", " | \"newton.linear.tolerance\": 1e-5, # tolerance for linear solver\n", " | \"newton.linear.errormeasure\": \"absolute\", # or \"relative\" or \"residualreduction\"\n", " | \"newton.linear.preconditioning.method\": \"jacobi\", # (see table below)\n", " | \"newton.linear.preconditioning.hypre.method\": \"boomeramg\", # \"pilu-t\" \"parasails\"\n", " | \"newton.linear.preconditioning.iteration\": 3, # iterations for preconditioner\n", " | \"newton.linear.preconditioning.relaxation\": 1.0, # omega for SOR and ILU\n", " | \"newton.linear.maxiterations\":1000, # max number of linear iterations\n", " | \"newton.linear.verbose\": False, # toggle linear iteration output\n", " | \"newton.linear.preconditioning.level\": 0} # fill-in level for ILU preconditioning\n", " | -----------------------------------------------\n", " | | Precondition | (x = parallel | s = serial) |\n", " | | method | numpy | istl | petsc |\n", " | |---------------|---------|---------|---------|\n", " | | amg-ilu | --- | x | --- |\n", " | | amg-jacobi | --- | x | --- |\n", " | | gauss-seidel | x | x | x |\n", " | | hypre | --- | --- | x |\n", " | | icc | --- | --- | x |\n", " | | ildl | --- | x | --- |\n", " | | ilu | --- | s | s |\n", " | | jacobi | x | x | x |\n", " | | kspoptions | --- | --- | x |\n", " | | lu | --- | --- | s |\n", " | | ml | --- | --- | x |\n", " | | none | x | x | x |\n", " | | oas | --- | --- | x |\n", " | | pcgamg | --- | --- | x |\n", " | | sor | x | x | x |\n", " | | ssor | x | x | x |\n", " | -----------------------------------------------\n", " | \n", " | The functionality of some of the preconditioners listed for petsc will\n", " | depend on the petsc installation.\n", " | \n", " | Method resolution order:\n", " | Scheme\n", " | pybind11_builtins.pybind11_object\n", " | builtins.object\n", " | \n", " | Methods defined here:\n", " | \n", " | __call__(...)\n", " | \n", " | __init__(...)\n", " | \n", " | dirichletIndices = _opDirichletIndices(self, id=None)\n", " | \n", " | inverseLinearOperator(...)\n", " | \n", " | jacobian(...)\n", " | \n", " | linear = _schemeLinear(self, assemble=True, parameters=None)\n", " | \n", " | setErrorMeasure(...)\n", " | \n", " | setQuadratureOrders(...)\n", " | \n", " | solve(scheme, target, rhs=None)\n", " | \n", " | ----------------------------------------------------------------------\n", " | Readonly properties defined here:\n", " | \n", " | cppIncludes\n", " | \n", " | cppTypeName\n", " | \n", " | dimRange\n", " | \n", " | domainSpace\n", " | \n", " | parameterHelp\n", " | \n", " | rangeSpace\n", " | \n", " | space\n", " | \n", " | ----------------------------------------------------------------------\n", " | Data descriptors defined here:\n", " | \n", " | __dict__\n", " | \n", " | ----------------------------------------------------------------------\n", " | Data and other attributes defined here:\n", " | \n", " | LinearInverseOperator =