Loading files
This file is pinned to the selected release. Download the original to open it in your design tool.
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(complementary:preconditioners)=\n",
"# Preconditioners in NGSolve\n",
"\n",
"Preconditioners are approximate inverses which are used within iterative methods to solve linear or non-linear equations.\n",
"\n",
"Here are some built-in preconditioners in NGSolve:\n",
"\n",
"* Jacobi (`local`) and block Jacobi \n",
"* Direct solvers using sparse factorization (`direct`)\n",
"* Geometric multigrid with different block-smoothers (`multigrid`)\n",
"* Algebraic multigrid preconditioner (`h1amg`)\n",
"* p-version element-level BDDC (`bddc`)\n",
"\n",
"This tutorial quickly shows how to use some of these within a solver."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from ngsolve import *\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A simple test problem \n",
"\n",
"In order to experiment with various preconditioners, \n",
"let us define a simple Poisson solver with the name\n",
"of a preconditioner as the keyword argument `precond`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def SolveProblem(h=0.5, p=1, levels=1, \n",
" condense=False,\n",
" precond=preconditioners.Local, preflags={}):\n",
" \"\"\"\n",
" Solve Poisson problem on l refinement levels.\n",
" PARAMETERS:\n",
" h: coarse mesh size\n",
" p: polynomial degree \n",
" l: number of refinement levels\n",
" condense: if true, perform static condensations\n",
" precond: preconditioner class (formerly: name of a built-in preconditioner)\n",
" Returns: (ndof, niterations)\n",
" List of tuples of number of degrees of freedom and iterations\n",
" \"\"\"\n",
" mesh = Mesh(unit_square.GenerateMesh(maxh=h))\n",
" # mesh = Mesh(unit_cube.GenerateMesh(maxh=h))\n",
" fes = H1(mesh, order=p, dirichlet=\"bottom|left\")\n",
" \n",
" u, v = fes.TnT() \n",
" a = BilinearForm(grad(u)*grad(v)*dx, condense=condense)\n",
" f = LinearForm(v*dx)\n",
" gfu = GridFunction(fes)\n",
" Draw(gfu)\n",
" c = precond(a,**preflags)\n",
"\n",
" steps = []\n",
" \n",
" for l in range(levels):\n",
" if l > 0: \n",
" mesh.Refine()\n",
" fes.Update()\n",
" gfu.Update()\n",
"\n",
" with TaskManager():\n",
" a.Assemble()\n",
" f.Assemble()\n",
"\n",
" # Conjugate gradient solver\n",
" inv = CGSolver(a.mat, c.mat, maxsteps=1000)\n",
"\n",
" # Solve steps depend on condense \n",
" if condense:\n",
" f.vec.data += a.harmonic_extension_trans * f.vec\n",
" \n",
" gfu.vec.data = inv * f.vec\n",
" \n",
" if condense:\n",
" gfu.vec.data += a.harmonic_extension * gfu.vec\n",
" gfu.vec.data += a.inner_solve * f.vec\n",
" steps.append ( (fes.ndof, inv.GetSteps()) )\n",
" if fes.ndof < 15000:\n",
" Redraw()\n",
" return steps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `Preconditioner` registers itself to the `BilinearForm`. Whenever the `BilinearForm` is re-assembled, the `Preconditioner` is updated as well."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The `local` preconditioner \n",
"\n",
"The `local` preconditioner is a simple Jacobi preconditioner. \n",
"The number of CG-iterations with the local preconditioner is proportional to $h^{-1} \\sim 2^l$:\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[(8, 4)]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"SolveProblem(precond=preconditioners.Local)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[(8, 4),\n",
" (21, 6),\n",
" (65, 16),\n",
" (225, 31),\n",
" (833, 61),\n",
" (3201, 119),\n",
" (12545, 236),\n",
" (49665, 452),\n",
" (197633, 906)]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res_local = SolveProblem(levels=9, precond=preconditioners.Local, preflags={\"GS\":True})\n",
"res_local"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multigrid preconditioner \n",
"\n",
"A geometric multigrid `Preconditioner` uses the sequence of refined meshes to define a preconditioner yielding optimal iteration numbers (and complexity). It uses a direct solve on the coarsest level, and block Gauss-Seidel smoothers on the refined levels."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[(8, 2),\n",
" (21, 4),\n",
" (65, 7),\n",
" (225, 8),\n",
" (833, 8),\n",
" (3201, 8),\n",
" (12545, 8),\n",
" (49665, 8),\n",
" (197633, 8)]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# res_mg = SolveProblem(levels=9, precond=\"multigrid\")\n",
"res_mg = SolveProblem(levels=9, precond=preconditioners.MultiGrid)\n",
"res_mg"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAGhCAYAAABVk3+7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABHB0lEQVR4nO3dd3wUdeLG8c9uegcSSCGhKS2kKSRY8ABFEREERFHvFLHdeVg59UDsiujpIZaoP70Ty+mJhWb34EAUkWrovbcEAukhbXd+fyxEOAgkYbOz5Xm/XnklO7vZfYZJsg/znfmOxTAMAxERERE3YDU7gIiIiMhRKiYiIiLiNlRMRERExG2omIiIiIjbUDERERERt6FiIiIiIm5DxURERETchr/ZARrKbrezd+9eIiIisFgsZscRERGRejAMg5KSEhISErBa694v4nHFZO/evSQlJZkdQ0RERBph165dJCYm1nm/xxWTiIgIwLFikZGRJqcRERGR+iguLiYpKan2fbwuHldMjg7fREZGqpiIiIh4mNMdhqGDX0VERMRtqJiIiIiI21AxEREREbfhcceY1JfNZqO6utrsGB4pMDDwlKdyiYiINBWvKyaGYZCbm0thYaHZUTyW1Wqlffv2BAYGmh1FRER8jNcVk6OlpFWrVoSGhmoStgY6OoHdvn37aNOmjf79RETEpbyqmNhsttpSEh0dbXYcj9WyZUv27t1LTU0NAQEBZscREREf4lUHEhw9piQ0NNTkJJ7t6BCOzWYzOYmIiPgaryomR2n44czo309ERMzi8mJSWFhIjx49yMjIICUlhbffftvVEURERMRNufwYk4iICObPn09oaChlZWWkpKQwbNgwHRMiIiJispW7C5n49XrGXdGFtMRmpmRw+R4TPz+/2mNAKisrMQwDwzBcHcPt9OnTh/vuu88lr7V9+3YsFgs5OTkueT0REfEM05bvYeHWg0xbvse0DA0uJvPnz2fQoEEkJCRgsViYMWPGCY/Jzs6mXbt2BAcH07NnTxYvXnzc/YWFhaSnp5OYmMiDDz5ITExMo1egqazcXcj1b/3Cyt2FZkcRERFpMrsLylm1u4jVe4qYtWIvALNy9rB6TxGrdhexu6DcpXkaPJRTVlZGeno6t9xyC8OGDTvh/qlTpzJmzBjefPNNevbsyeTJk+nfvz8bNmygVatWADRr1owVK1aQl5fHsGHDGD58OLGxsSd9vcrKSiorK2tvFxcXNzRyoxzbGs3anSUiIuJsNTY7uwsOszW/lK0Hynjmq3UnPOZQeTVXvvpT7e3tzw10Wb4GF5MBAwYwYMCAOu+fNGkSt99+O6NGjQLgzTff5KuvvuKdd95h7Nixxz02NjaW9PR0fvzxR4YPH37S55s4cSJPPvlkQ2MCjllgD1fX/5TXPYWHKSyvwoLlt9a4Yi9XpsVjYNAsNJDWzULq9VwhAX6NPruloKCAe++9ly+++ILKykp69+7NK6+8QseOHWsfs2DBAsaPH8/ixYsJCgoiKyuLjz/+mObNm/Ptt9/yzDPPsHr1avz8/Dj//PN5+eWXOeussxqVR0REPM+hsiq2HnCUj635ZY6v88vYcbCMalv9DqHwt1p48Zr0Jk76P6/pzCerqqpi2bJljBs3rnaZ1WqlX79+LFy4EIC8vDxCQ0OJiIigqKiI+fPnc+edd9b5nOPGjWPMmDG1t4uLi0lKSqpXnsPVNpIf+66Ra+NwqKyK4W8ubPD3rX2qP6GBjfvnvfnmm9m0aROzZs0iMjKSv/71r1xxxRWsXbuWgIAAcnJyuOSSS7jlllt4+eWX8ff3Z+7cubXzjpSVlTFmzBjS0tIoLS3lscceY+jQoeTk5OgaOCIiXqSyxsbOg+VsOVBWuwfkaAEpLK/7enFB/lbax4TRoWUYHWLCCfS3MOk/m0543IzRF5LSOqopV+EETi0m+fn52Gy2E4ZlYmNjWb9+PQA7duzgjjvuqD3o9e677yY1NbXO5wwKCiIoKMiZMd3a0UKyYMECLrjgAgA+/PBDkpKSmDFjBtdccw1/+9vf6NGjB6+//nrt93Xr1q3266uvvvq453znnXdo2bIla9euJSUlxTUrIiIiTmEYBvtLKtlydO/HMSVkd0E59lPs/EiICqZDy/AjBSSs9uuEqBCs1t/26q/eU8Sk/2zCYgHDoPazGVx+unBWVpbLzgYJCfBj7VP9G/Q9a/cWn3QPyWd/Op/khMgGvXZjrFu3Dn9/f3r27Fm7LDo6ms6dO7NunWMcMCcnh2uuuabO59i0aROPPfYYixYtIj8/H7vdDsDOnTtVTERE3FR5Vc3xwy4HytiW7/gorayp8/vCg/xPKB4dYsJpHxNGSGD93ouiwwNpGR5EfLNgRmQmMXXJLvYVVhAd7vqLuTq1mMTExODn50deXt5xy/Py8oiLi3PmS9WLxWJp8HBK8JFC8b+tMTjAr9FDM84WEnLq41wGDRpE27Ztefvtt0lISMBut5OSkkJVVZWLEoqIyMnY7AZ7Cw+z5UAp2/KP3/uxr6iizu/zs1pIah7iKB5HCkj7mDDOahlGy4igM56xOz4qhJ/G9iXQz4rFYuGGrDZU2ewE+TfuP9lnwqnvtIGBgXTv3p05c+YwZMgQwHG12jlz5nDXXXc586WajNmtsWvXrtTU1LBo0aLaoZyDBw+yYcMGkpOTAUhLS2POnDknPSj46GPffvttLrroIgB++umnEx4nIiKndiaTjRWVV/92zEf+b0Mw2w6WUVVjr/P7WoQFHikeYbSPcez9OKtlGG1ahBHo37THCB5bQiwWiymlBBpRTEpLS9m8eXPt7W3btpGTk0OLFi1o06YNY8aMYeTIkfTo0YOsrCwmT55MWVlZ7Vk6jZWdnU12dnaTX1jO7NbYsWNHrrrqKm6//Xb+7//+j4iICMaOHUvr1q256qqrAMcBwampqfz5z3/mT3/6E4GBgcydO5drrrmGFi1aEB0dzVtvvUV8fDw7d+484WwoERE5vdNNG1Fts7PzUPlvB5weU0IOltW9hzrQz0rb6FDHkMsxe0DOahlGs1DXD524mwYXk6VLl9K3b9/a20fPmBk5ciTvvvsuI0aM4MCBAzz22GPk5uaSkZHBt99+W+c8JfU1evRoRo8eTXFxMVFRTXuEsNmtccqUKdx7771ceeWVVFVV8bvf/Y6vv/6agIAAADp16sT333/Pww8/TFZWFiEhIfTs2ZPrr78eq9XKxx9/zD333ENKSgqdO3fmlVdeoU+fPi5dBxERT7S7oJyCsmosFvjiyLQRM3P20Dkugj0Fh8kvreBASRXb8svYeaicmlMceRobGUSHI3s9jh77cVZMOK2bh+Bn1cVS62IxPGw++KPFpKioiMjI4w9GraioYNu2bbRv357g4GCTEno+/TuKiK9qN/arBj0+NNDvyGm34bVDMB1iwmnfMozwIPc4LtFdnOr9+1j6VxMRETni2aGpjJ++ipP9j90CXNQxhkuTY2v3gMRFBp/xgadyPBUTERHxeYZhMCNnD5P+s+GkpQTgi7t7uXyyMV+kYiIiIj5ty4FSHp2xmp+3HAQgsXkIuwsOu8VkY77IY4qJq87KERER31BRbeP1eVt4c96WI2dfWrnnko4MSovn6jcWusVkY75IB7/KCfTvKCLebv7GAzw6czU7DpYD0KdzS54anEKb6FDAcQ2ao9NGGIZh2mRj3kQHv4qIiPyP/cUVPPXlWr5cuQ9wnNL7xKBuXJ4Sd9xBrGZPG+HLVExERMTr2ewG//plBy9+t4GSyhqsFrj5gvaMuayTTut1M9oaIiLi1VbtLmL8jFWs3F0EQHpiFBOGpuoMGzelYiIiIl6puKKaSd9v5P2F27EbEBHsz0OXd+GGrDaaedWNqZiIiIhXMQyDr1bt46kv1rK/pBKAqzISGD+wK60idEC/u/OYYuLy04VLcmHpFOgxCiLiXPOaIiJyRnYcLOPRmWuYv/EAAO1jwnj6qhR6dYwxOZnUV9NeQ9mJRo8ezdq1a1myZIlrXrAkF354zvHZBfr06cPdd9/NfffdR/PmzYmNjeXtt9+uvTJzREQEZ599Nt98803t98yaNYuOHTsSHBxM3759ee+997BYLBQWFroks4iIu6issfHqnE1c9tJ85m88QKCflfv6deSbey9SKfEwHrPHpFEMA6rLG/e9NYd/+1xV1vDvDwh1TBfYAO+99x4PPfQQixcvZurUqdx5551Mnz6doUOH8vDDD/PSSy9x4403snPnTvLy8hg+fDj33nsvt912G7/++isPPPBAw3OKiHi4hVsO8siMVWw54Phb3evsGJ4ekkL7mDCTk0ljePcEa1Vl8GyCOUEf3guB9f+l6NOnDzabjR9//BEAm81GVFQUw4YN4/333wcgNzeX+Ph4Fi5cyIwZM/jqq69YtWpV7XM88sgjTJgwgYKCApo1a9bo6JpgTUQ8QX5pJc9+tY5pv+4BICY8iEev7Mrg9ARdWM8NaYI1D5SWllb7tZ+fH9HR0aSmptYui42NBWD//v1s2LCBzMzM474/KyvLNUFFRExktxt8vGQXz3+7nqLD1Vgs8IeebXmgf2eiQgLMjidnyLuLSUCoY89FfZXmOT4AclfB1w/CFS9A3JFyEB7r+KjvazdQQMDxv1AWi+W4ZUf/B2C32xv83CIi3mDt3mIembGK5TsLAeiWEMmEoalkJDUzNZc4j3cXE4ulQcMptOjg+ADwD3F8TsyChAynRztTnTt35uuvvz5umcsODBYRcbGyyhomz97IOwu2Y7MbhAf5M+bSTtx0flv8/TzmPA6pB+8uJl7sj3/8I5MmTeKvf/0rt956Kzk5Obz77rsAGlsVEa/y3Zpcnpi1hn1FFQAMTI3n0SuTiYvSMXDeSDWzLhFx0Hus285h0r59ez777DOmTZtGWloab7zxBuPHjwcgKCjI5HQiImdud0E5t723hD9+sIx9RRUktQhhyqhMsn9/rkqJF/OYPSYun2AtIg76jnPNawHz5s07Ydn27dtPWHbsSVSDBw9m8ODBtbcnTJhAYmKizqQREY9WbbPzz5+28fLsTRyuthHgZ+GO33Xgrr4dCQnUVX69nccUk9GjRzN69Oja040EXn/9dTIzM4mOjmbBggW88MIL3HXXXWbHEhFptCXbD/HI9NVsyCsBIKt9CyYMSaFjbITJycRVPKaYyIk2bdrEM888w6FDh2jTpg1/+ctfGDfOdXt5REScpaCsiue+Wc/UpbsAaBEWyMNXdOXqc1vruDkfo2LiwV566SVeeukls2OIiDSaYRh8tmw3z369joLyagCuz0riof5daB4WaHI6MYOKiYiImGJTXgnjZ6xm8bZDAHSJi+CZISn0aNfC5GRiJq8sJh42y77b0b+fiDSlw1U2Xv3vJt6av5Uau0FIgB/3X9qRURe2J0Bzkvg8ryomR2dJLS8vJyQkxOQ0nquqqgpwTIsvIuJMc9fv59GZq9ld4LhQ6qXJsTwxuButm+lvtjh4VTHx8/OjWbNm7N+/H4DQ0FAdNNVAdrudAwcOEBoair+/V/14iIiJ9hUd5qkv1vLN6lwAEqKCeWJwNy7r5p5zRYl5vO6dJy7O8UN+tJxIw1mtVtq0aaNSJyJnrMZm572FO5j0/QbKqmz4WS3c1qs991zSkbAgr3sLEifwmJ+K+k6wZrFYiI+Pp1WrVlRXV7sonXcJDAzEatU4r4icmZxdhTw8bRVr9xUD0L1tc54ZkkLX+LoveS9iMTzsSMejE6wVFRURGakfbhERd1N0uJoXvlvPh4t2YhgQFRLAuAFduLZHElar9sT6qvq+f3vMHhMREXFvhmEwa8Venv5yHfmllQBcfW4iD1/RhehwXcNL6kfFREREztjWA6U8OnM1CzYfBOCslmE8MySV88+KNjmZeBoVExERabSKahtvzNvCG/O2UGWzE+Rv5Z5LOnL7RR0I9NexatJwKiYiItIoP246wKMzVrP9YDkAvTu15OmrUmgTHWpyMvFkKiYiInJKK3cXMvHr9Yy7ogtpic3YX1LBM1+uY9aKvQDERgbx+KBuDEiJ0zQDcsZUTERE5JSmLd/Dwq0H+XzZHlbsKuRv322gpKIGqwVGXtCOMZd2IiI4wOyY4iVUTERE5AS7C8opKKvGYoEvjuwZ+XDRDmrsjhkmusZF8MI16aS0jjIzpnghFRMRETlBr+fnnrDsaCkBWJdbolIiTUKHTIuIyAleujaduuZC87damDwiw6V5xHdoj4mIiBxn16FyZq3Yi72OecFnjL5Qe0ukyXhMManvtXJERKRxqm12/vHjNl6es5GKajv+Vgs1dgOLBQyD2s8iTUnXyhEREZZuP8TD01exMa8UgPM6tODeSzpxz79/Jb5ZMCMyk5i6ZBf7CiuYdfeFxEeFmJxYPE19379VTEREfFhheRXPfbOej5fsAqBFWCCPDOzK0HNaY7FYqKyxEehnxWKxYBjGkdld/UxOLZ5IF/ETEZE6GYbBtOV7mPD1Og6VVQFwXWYSYwd0oVloYO3jji0hFotFpUSanIqJiIiP2by/lEdmrOKXrYcA6BwbwYShKfRo18LkZCIqJiIiPqOi2sbrczfzxg9bqLYZBAdYufeSTtx2UXsC/DR7hLgHFRMRER/w46YDPDJjNTuOXHCvb+eWPHVVCkktdME9cS8qJiIiXuxkF9x7YlA3LtcF98RNqZiIiHghu93gw8U7+du364+74N5fLutMeJD+9Iv70k+niIiXWbO3iPHTV5OzqxCAtMQonh2aqtlaxSOomIiIeImyyhpe+s9Gpvy8HZvdIDzInwf7d+YP57XFr64L34i4GRUTEREv8N2aXJ6YtYZ9RRUADEyL57Erk4mNDDY5mUjDqJiIiHiw3QXlPDFrLbPX5QGQ1CKEp69KoU/nViYnE2kcFRMREQ9UbbMzZcE2XvrPJg5X2wjws3DH7zpwV9+OhARqdlbxXComIiIeZtmOAsZPX8X63BIAstq1YMLQFDrGRpicTOTMqZiIiHiIovJqnv9uPf9evBPDgOahAYy7oivXdE/UnCTiNTymmGRnZ5OdnY3NZjM7ioiISxmGwcycvTzz1VrySx0X3LumeyLjruhKi7DA03y3iGexGIZhmB2iIep72WQREW+w9UApj85czYLNBwE4u1U4E4ak0LNDtMnJRBqmvu/fHrPHRETEl1RU23jzhy28PncLVTY7Qf5W7rmkI7df1IFAf11wT7yXiomIiJv5eXM+j8xYzdb8MgB+16klT1/VjbbRYSYnE2l6KiYiIm4iv7SSCV+tY/qvewBoGRHE44OSGZgar4NbxWeomIiImMxuN/h4yS6e+2YdxRU1WCxw03lt+Uv/zkQGB5gdT8SlVExEREy0bl8x46evYvnOQgC6JUTy7NBU0pOamZpLxCwqJiIiJiivquHl2Zv4x0/bsNkNwgL9+Mtlnbnp/Lb4++ngVvFdKiYiIi42e20ej89aw57CwwAMSInj8UHdiIvSBfdEVExERFxkb+FhnvxiDd+tcVxwr3WzEJ4e0o2Lu8SanEzEfaiYiIg0sRqbnXd/3s5L/9lIWZUNf6uF2y7qwD2XnE1ooP4MixxLvxEiIk0oZ1chD09bxdp9xQB0b9ucCUNT6BKnmatFTkbFRESkCRQdrubF7zbwr0U7MAyICglg3IAuXNsjCatVc5KI1EXFRETEiQzD4IuV+3j6y7UcKKkEYNg5rXl4YFdiwoNMTifi/lRMREScZHt+GY/OXM2Pm/IB6NAyjGeGpHDBWTEmJxPxHComIiJnqLLGxls/bOXVuZupqrET6G/lrr5n88feHQjy9zM7nohHUTEREWmglbsLmfj1esZd0YWyShuPzFjFlgOOC+71OjuGp4ek0D5GF9wTaQwVExGRBpq2fA8Ltx7k/qkr2HKgFICY8CAevbIrg9MTdME9kTOgYiIiUg+7C8opKKvGMAw+XbYLoLaUDEiJ455LOtI1XqcAi5wpFRMRkXro9fzcOu/7ZnUu36zOZftzA12YSMQ76UpRIiKnUVlj4/JucXXe72+1MHlEhusCiXgx7TERETmFJdsPMW7aKjbvL63zMTNGX0hK6ygXphLxXh6zxyQ7O5vk5GQyMzPNjiIiPqDocDUPT1/FNW8uZPP+UmLCA3mof2cAjh7bqmNcRZzPYhiGYXaIhiguLiYqKoqioiIiI3WgmYg4l2EYfLM6l8dnramdufW6zCTGDujC4Wobg19dQHyzYEZkJjF1yS72FVYw6+4LiY8KMTm5iHur7/u3iomIyBF7Cw/z2MzVzF63H4AOMWE8OyyV8zpE1z6mssZGoJ8Vi8WCYRhU2eyaRE2kHur7/q1jTETE59nsBu8v3M6L322grMpGgJ+FO/uczZ/7nEVwwPGl49gSYrFYVEpEnEzFRER82tq9xYybtpIVu4sA6NG2OROHpdIxNsLkZCK+ScVERHxSRbWNybM38faPW7HZDSKC/PnrgC7ckNUGq1VHtYqYRcVERHzOT5vyGT9jFTsOlgOOmVufGNyN2Mhgk5OJiIqJiPiMQ2VVPPPlWqb9ugeAuMhgnh6SwqXJsSYnE5GjVExExOsZhsH0X/fw9JdrKSivxmKBkee34y+XdSIiOMDseCJyDBUTEfFqOw6WMX76an7anA9Al7gIJg5L5Zw2zU1OJiIno2IiIl6p2mbnHz9uY/LsjVTW2Anyt3Jvv47cflEHAvw8ZtJrEZ+jYiIiXidnVyFjP1/J+twSAC48O5oJQ1JpFxNmcjIROR0VExHxGqWVNbz43QbeW7gdw4DmoQE8MjCZYee2xqIL24h4BBUTEfEKs9fm8ejM1ewrqgBg2DmtGT+wK9HhQSYnE5GGUDEREY+2v7iCJ75Yw9ercgFIahHChCGp/K5TS5OTiUhjqJiIiEey2w3+vWQnz32znpKKGvysFm67qD33XdKJkEBdv0bEU6mYiIjH2by/hHHTVrFkewEAaYlRTByWSreEKJOTiciZUjEREY9RWWPj9blbeH3eZqptBqGBfjxwWWdGXtAOP13fRsQrqJiIiEdYtPUg46avYuuBMgAu7tKKp4ek0LpZiMnJRMSZVExExK0VlVfz3Lfr+PfiXQDEhAfx5OBuXJEap1OARbyQiomIuCXDMPhq1T6emLWW/NJKAK7PasPYy7sQFarr24h4KxUTEXE7ewoP8+iM1fx3/X4AzmoZxsRhaWS1b2FyMhFpaiomIuI2bHaDd3/ezt+/30B5lY0APwt/7nM2f+57FkH+OgVYxBeomIiIW1izt4hx01axcncRAJntmjNxWCpnt4owOZmIuJKKiYiY6nCVjclzNvKPH7dhsxtEBPszbkBXrstMwqpTgEV8joqJiJhm/sYDjJ+xil2HDgMwMDWexwcl0yoy2ORkImIWFRMRcbmDpZU889U6pv+6B4CEqGCeuiqFfsmxJicTEbOpmIiIyxiGwefL9/DMV2spLK/GYoGbL2jHXy7rTHiQ/hyJiIqJiLjI9vwyHp6+ip+3HASga3wkE4elkpHUzNxgIuJWVExEpElV2+y8NX8rr8zZRGWNnSB/K/df2olbe7UnwM9qdjwRcTMqJiLSZH7dWcC4aatYn1sCQK+zY5gwNIW20WEmJxMRd6ViIiJOV1JRzYvfbeD9X3ZgGNAiLJBHBnZl6DmtdX0bETkll+9H3bVrF3369CE5OZm0tDQ+/fRTV0cQESdbubuQ69/6hZW7C/l+TS6XTprPewsdpWTYua2ZPaY3w85NVCkRkdNy+R4Tf39/Jk+eTEZGBrm5uXTv3p0rrriCsDDt2hXxVNOW72Hh1oPc/dGv7DhUDkCbFqE8OzSVXh1jTE4nIp7E5cUkPj6e+Ph4AOLi4oiJieHQoUMqJiIeZndBOQVl1RiGwadLdwGw41A5VgsMOzeRP/U+i7NbhZucUkQ8TYOHcubPn8+gQYNISEjAYrEwY8aMEx6TnZ1Nu3btCA4OpmfPnixevPikz7Vs2TJsNhtJSUkNDi4i5ur1/FwGvfYTg7MXUFZlq11uN+CzZbvpN+kHE9OJiKdqcDEpKysjPT2d7Ozsk94/depUxowZw+OPP87y5ctJT0+nf//+7N+//7jHHTp0iJtuuom33nrrlK9XWVlJcXHxcR8iYq4am50r0+LrvN/famHyiAzXBRIRr2ExDMNo9DdbLEyfPp0hQ4bULuvZsyeZmZm89tprANjtdpKSkrj77rsZO3Ys4Cgbl156Kbfffjs33njjKV/jiSee4MknnzxheVFREZGRkY2NLiKNtHZvMQ99voLVe+r+T8KXd/cipXWUC1OJiLsrLi4mKirqtO/fTj0rp6qqimXLltGvX7/fXsBqpV+/fixcuBBwTEl98803c/HFF5+2lACMGzeOoqKi2o9du3Y5M7KI1FNljY2/f7+Bwa/9xOo9xUQG+3PvJWcDcPRkG510IyJnyqkHv+bn52Oz2YiNPf5CXLGxsaxfvx6ABQsWMHXqVNLS0mqPT/nggw9ITU096XMGBQURFBTkzJgi0kDLdxbw0Gcr2by/FIDLu8Xx1JBu2OwGHy3aRXyzYEZkJjF1yS72FVYQHR5ocmIR8VQuPyunV69e2O12V7+siDRCeVUNL363kSk/b8MwICY8kKeuSuGK1N+OL/lpbF8C/axYLBZuyGpDlc1OkL+fialFxJM5tZjExMTg5+dHXl7eccvz8vKIi4tz5kuJSBP7eXM+Y6etYueReUmGnduaRwcm0zzs+L0hx5YQi8WiUiIiZ8Spx5gEBgbSvXt35syZU7vMbrczZ84czj//fGe+lIg0keKKasZ+vpIb/rGInYfKSYgKZsqoTCZdm3FCKRERcbYG7zEpLS1l8+bNtbe3bdtGTk4OLVq0oE2bNowZM4aRI0fSo0cPsrKymDx5MmVlZYwaNeqMgmZnZ5OdnY3NZjv9g0WkUWavzWP8jFXkFVcCcON5bfnrgC6EB+myWiLiGg0+XXjevHn07dv3hOUjR47k3XffBeC1117jhRdeIDc3l4yMDF555RV69uzplMD1Pd1IROrvYGklT3yxli9W7AWgfUwYzw1LpWeHaJOTiYi3qO/79xnNY2IGFRMR5zEMg1kr9vLkF2s5VFaF1QK3/64D9/frRHCAjhUREeep7/u39s+K+KjcogrGT1/FnPWOWZm7xEXwt+FppCU2MzeYiPg0FRMRH2MYBh8v2cWzX62jpLKGAD8Ld/XtyJ19ziLQ36nHw4uINJiKiYgP2XGwjLGfr2Lh1oMAZCQ142/D0+gUG2FyMhERB48pJjorR6TxbHaDKQu28eL3G6iothMcYOWByzoz6sL2+Fk1j7yIuA8d/Cri5TbmlfDQZyvJ2VUIwPkdonnu6lTaRoeZG0xEfIoOfhXxcVU1dt78YQuv/ncT1TaDiCB/Hh7Ylesyk7Doansi4qZUTES80MrdhTz02UrW55YAcEmXVjwzNIX4qBCTk4mInJqKiYgXqai28dLsjbw9fyt2A5qHBvDE4G4MTk/QXhIR8QgqJiJeYvG2Q/z185Vsyy8DYFB6Ao8PSiYmPMjkZCIi9adiIuLhSitreP6b9Xzwyw4AWkUEMWFoKpcmx5qcTESk4VRMRDzYvA37GT99NXsKDwNwXWYS467oSlRIgMnJREQax2OKieYxEflNYXkVT325lmnL9wCQ1CKE54alceHZMSYnExE5M5rHRMTDfL1qH4/NXE1+aRUWC4y6oD0P9O9EaKDH/D9DRHyQ5jER8TL7Syp4bMYavl2TC8DZrcJ5/uo0urdtbnIyERHnUTERcXOGYfD58j08/eVaig5X42+1cGefs7jr4rMJ8vczO56IiFOpmIi4sd0F5Tw8fTXzNx4AIKV1JM9fnUa3hCiTk4mINA0VExE3ZLcb/GvRDp7/Zj1lVTYC/a3c168jd1zUAX8/q9nxRESajIqJiJvZcqCUsZ+vZMn2AgAy2zXnuavTOKtluMnJRESanoqJiJuosdl5+8dtvDR7I1U1dkID/Rg7oAt/6NkWq1XTyYuIb/CYYqJ5TMSbrd1bzEOfr2D1nmIALuoYw8RhqSQ2DzU5mYiIa2keExETVdbYeO2/m3lj3hZq7AaRwf48emUyw7sn6qJ7IuJVNI+JiJtbvrOAhz5byeb9pQBc3i2Op4Z0o1VEsMnJRETMo2Ii4mLlVTW8+N1Gpvy8DcOAmPBAnroqhStS482OJiJiOhUTERdasDmfsdNWsuuQ46J7w85tzaMDk2keFmhyMhER96BiItJEVu4uZOLX6xl3RRfaRocx8et1fLxkFwAJUcFMGJZK386tTE4pIuJeVExEmsi05XtYuPUgL8/exOq9ReQVVwJw43lteejyzkQEB5icUETE/aiYiDjR7oJyCsqqsVhgZs4eAOas3w9AQrNgHh7QlSvTE8yMKCLi1lRMRJyo1/Nz67xvb2EFd/37VxUTEZFT0EU3RJzo+atTqWv2EX+rhckjMlwZR0TE43jMHhPN/CrubtmOArLnbqGuGQtnjL6QlNa6KrCIyKl4zB6T0aNHs3btWpYsWWJ2FJHjVNvs/P37DVzz5s/sPFROTLjj1N+jE7dqAlcRkfrzmGIi4o427y9h2Os/8+p/N2M3YNg5rfno9vNoGR5EausoJgxNIbV1FC3Dg4gO11wlIiKno2vliDSC3W7w/sLtTPxmPZU1dqJCAnh2aCoD0xyzt1bW2Aj0s2KxWDAMgyqbnSB/P5NTi4iYR9fKEWkiuUUVPPjZCn7clA84rgT8wvB04qJ+u8bNsSXEYrGolIiI1JOKiUgDfLlyL+Onr6bocDVB/lbGD+zKjee11ZWARUScRMVEpB6KDlfz+MzVzMjZC0Bq6yheGpHB2a3CTU4mIuJdVExETuPnzfn85dMV7CuqwGqBu/qezd2XdCTAT8eOi4g4m4qJSB0qqm28+N0G/vHTNgDaRocy6doMurdtbnIyERHvpWIichJr9hZx/9QcNuaVAnB9VhseGdiVsCD9yoiINCX9lRU5hs1u8Nb8rUz6zwaqbQYx4YE8f3Ual3SNNTuaiIhPUDEROWLXoXL+8skKFm8/BMClybE8NyyV6PAgk5OJiPgOjykmulaONBXDMPhs2W6e/GItpZU1hAX68figblzTI1GnAYuIuJhmfhWfdqisinHTVvLdmjwAerRtzqRrM2gTHWpyMhER76KZX0VOY+76/Tz42UrySyvxt1q4/9JO/Kn3WfhZtZdERMQsKibic8qrapjw1To+XLQTgI6twnlpRAYpraNMTiYiIiom4lN+3VnAmE9WsC2/DIBbLmzPQ5d3JjhA17IREXEHKibiE6ptdl7772Zem7sZm90gLjKYF69Jp1fHGLOjiYjIMVRMxOttOVDKmKk5rNhdBMDg9ASeviqFqNAAk5OJiMj/UjERr2UYBv9atJMJX62lotpOZLA/Tw9J4aqM1mZHExGROqiYiFfaX1zBQ5+vZN6GAwBceHY0LwxPJ6FZiMnJRETkVFRMxOt8s2ofD09fRUF5NYH+VsZe3oWbL2iHVacBi4i4PRUT8RrFFdU8MWsN05bvASA5PpKXr8ugY2yEyclERKS+VEzEKyzaepAxn6xgT+FhrBb4U++zuK9fJwL9rWZHExGRBlAxEY9WWWNj0vcbeevHrRgGJLUIYdK1GWS2a2F2NBERaQQVE/FY63OLue/jHNbnlgBwbY9EHhvUjfAg/ViLiHgq/QUXj2O3G/zzp2288N0Gqmx2WoQFMnFYKv27xZkdTUREzpDHFJPs7Gyys7Ox2WxmRxET7Sk8zF8+yeGXrYcAuKRLK567Oo2WEUEmJxMREWewGIZhmB2iIep72WTxLoZhMCNnD4/NWENJZQ0hAX48emUy12clYbHoNGAREXdX3/dvj9ljIr6rsLyK8dNX89WqfQCc06YZL12bQbuYMJOTiYiIs6mYiFv7YeMBHvx0BftLKvG3Wrj3ko7c2ecs/P10GrCIiDdSMRG3dLjKxnPfrOO9hTsA6NAyjMkjMkhLbGZuMBERaVIqJuJ2Vu4u5P6pOWw5UAbAyPPbMnZAV0IC/UxOJiIiTU3FREy1cnchE79ez7grupAcH8kb87bw8pxN1NgNWkUE8cI16fTu1NLsmCIi4iIqJmKqacv3sHDrQd77eTtb88v4dWchAANT43lmSArNwwLNDSgiIi6lYiIut7ugnIKyaiwW+GLFXgA+P3LhvdAAP/7SvxO3XNhepwGLiPggFRNxuV7Pz63zvvJqG09/uY5be3VwYSIREXEXOudSXG7yiAz8rCffG+JvtTB5RIZrA4mIiNtQMRGXG5yewO86xpz0vhmjL2TIOa1dnEhERNyFiom4lGEYTPh6HXM3HADg6H4THU4iIiKgY0zExd78YSv//GkbABHB/rSPCWNEZhJTl+xiX2EF0eE6C0dExJepmIjLfLJkF89/ux6A8Vd05aYL2hLoZ8VisXBDVhuqbHaC/DWJmoiIL1MxEZf4z9o8xk5bCcAfe3fg9t8df9aNxWJRKRERER1jIk1v8bZD3PXRcuwGDO+eyNjLu5gdSURE3JSKiTSpdfuKufW9JVTW2LmkSyueG5aqidNERKROKibSZHYdKmfkO4spqaihR9vmvHbDufj76UdORETqpncJaRL5pZXc+M9F7C+ppHNsBP8cmamrA4uIyGmpmIjTlVRUc/OUxWw/WE7rZiG8f2sWUaEBZscSEREPoGIiTlVZY+OPHyxj9Z5iWoQF8sGtWcRGBpsdS0REPITHFJPs7GySk5PJzMw0O4rUwWY3uH9qDj9vOUhYoB/vjsqkQ8tws2OJiIgHsRiGYZgdoiGKi4uJioqiqKiIyMhIs+PIEYZh8OjM1fzrl50E+FmYcnMWveq4Ho6IiPie+r5/e8weE3FvL8/ZxL9+2YnFAi+NyFApERGRRlExkTP2wS87mDx7EwBPDu7GlWkJJicSERFPpWIiZ+Srlft4bOZqAO65pCM3nd/O3EAiIuLRVEyk0RZszuf+qTkYBtzQsw339+todiQREfFwKibSKKt2F3HH+0upstkZkBLH01elaKp5ERE5Yyom0mDb8su4ecpiyqpsnN8hmsnXZeBnVSkREZEzp2IiDZJXXMGN/1zEwbIquiVE8tZN3Qny11TzIiLiHComUm9Fh6sZ+c5idhccpm10KO+OyiIiWFPNi4iI86iYSL1UVNu47b0lrM8toWVEEB/c0pOWEUFmxxIRES+jYiKnVWOzc9dHv7JkewERQf68NyqLNtGhZscSEREvpGIip2QYBg9PX8XsdXkE+lv5x8geJCfoUgAiItI0VEzklP723QY+WbobqwVeu/4cenaINjuSiIh4MRUTqdM/ftzKG/O2ADBxWCqXdYszOZGIiHg7FRM5qem/7uaZr9YB8GD/zozIbGNyIhER8QUqJnKCuRv28+CnKwEYdWE7/tznLJMTiYiIr1AxkeMs31nAn/+1nBq7wVUZCTw6MFlTzYuIiMuomEitTXkl3PLuEg5X2/hdp5a8MDwdq6aaFxERF1IxEQD2Fh7mpncWU1heTUZSM978w7kE+uvHQ0REXEvvPEJBWRU3/nMR+4oqOKtlGFNuziQ00N/sWCIi4oNUTHxceVUNo95dwpYDZcRHBfP+rT1pHhZodiwREfFRKiY+rKrGzp/+tZycXYU0Cw3g/VuyaN0sxOxYIiLiw1RMfJTdbvDgZyuYv/EAIQF+/HNkJh1jI8yOJSIiPk7FxAcZhsHTX61lZs5e/K0WXv/DuXRv29zsWCIiIiomvuj1eVuYsmA7AC9ck0bfzq3MDSQiInKEiomP+XjxTl74bgMAjwzsytBzEk1OJCIi8hsVEx/y3ZpcHp6+CoA7+5zFbRd1MDmRiIjI8VRMfMSirQe5+9+/Yjfg2h6JPNS/s9mRRERETqBi4gPW7i3mtveWUlVjp1/XWJ4dmqrr34iIiFtSMfFyOw+WM3LKYkoqa8hq14LXbjgHfz9tdhERcU96h/JiB0oqufGdRRwoqaRLXARvj+xBcICf2bFERETqpGLipUoqqrl5ymJ2HCwnsXkI792SRVRIgNmxRERETsmUYjJ06FCaN2/O8OHDzXh5r1dRbeOO95exZm8x0WGBfHBrT2Ijg82OJSIiclqmFJN7772X999/34yX9no2u8H9U3NYuPUgYYF+vDsqi/YxYWbHEhERqRdTikmfPn2IiNB1WZzNMAwenbmab1bnEuhn5a2bepCaGGV2LBERkXprcDGZP38+gwYNIiEhAYvFwowZM054THZ2Nu3atSM4OJiePXuyePFiZ2SV03hp9iY+WrQTiwVeGpHBhWfHmB1JRESkQRpcTMrKykhPTyc7O/uk90+dOpUxY8bw+OOPs3z5ctLT0+nfvz/79+9vVMDKykqKi4uP+5ATvb9wO6/M2QTAU1elMDAt3uREIiIiDdfgYjJgwACeeeYZhg4detL7J02axO23386oUaNITk7mzTffJDQ0lHfeeadRASdOnEhUVFTtR1JSUqOex5t9uXIvj89aA8B9/Tpy43ltTU4kIiLSOE49xqSqqoply5bRr1+/317AaqVfv34sXLiwUc85btw4ioqKaj927drlrLhe4adN+dw/NQfDgBvPa8u9l3Q0O5KIiEij+TvzyfLz87HZbMTGxh63PDY2lvXr19fe7tevHytWrKCsrIzExEQ+/fRTzj///JM+Z1BQEEFBQc6M6fFW7i5k4tfrGd69NY/OXEO1zWBgajxPDO6mqeZFRMSjObWY1Nfs2bPNeFmvMW35HhZuPcjynQVU1ti58OxoJo1Ix8+qUiIiIp7NqcUkJiYGPz8/8vLyjluel5dHXFycM1/K5+wuKKegrBqLBWbm7AGgssbO2S3DueeSjhwoqSSxeajJKUVERM6MU4tJYGAg3bt3Z86cOQwZMgQAu93OnDlzuOuuu5z5Uj6n1/NzT7p884FSRvzfLwBsf26gKyOJiIg4XYMPfi0tLSUnJ4ecnBwAtm3bRk5ODjt37gRgzJgxvP3227z33nusW7eOO++8k7KyMkaNGnVGQbOzs0lOTiYzM/OMnsdTTR6RUedQjb/VwuQRGa4NJCIi0gQshmEYDfmGefPm0bdv3xOWjxw5knfffReA1157jRdeeIHc3FwyMjJ45ZVX6Nmzp1MCFxcXExUVRVFREZGRkU55Tk9wqKyKYa8vYPvB8hPu+/LuXqS01gyvIiLivur7/t3gPSZ9+vTBMIwTPo6WEoC77rqLHTt2UFlZyaJFi5xWSnxVYXkVf/jHotpScnS/iU7AcXMluTB3ouOzt/OldQWtr6/wxfV2g3U25Vo5Un9Fh6u58Z+LWbuvmOahATQPDSA1MYoJQ1NIbR1Fy/AgosMDzY4pJ1OSCz885xt/1HxpXUHr6yt8cb3dYJ1NOV1Y6qe4opqb3lnMqj1FtAgL5OM7zqNtdCiBflYsFgs3ZLWhymYnyN/P7KhyMjVVjs9l+VC8z9wsTa0s/7fP3r6uoPX1Fb643kfXueawaREafIyJ2XzlGJPSyhpu+ucilu8spFloAB/ddh7JCd67vh7JboOSfVCwAwq2Q+EO2L8eDm11LC/PNzuhiEjjpF0H593p+DoizvFxhur7/u0xe0yys7PJzs7GZrOZHaXJlVXWMGrKYpbvLCQy2J9/3dpTpcQMhgGHCxyFo2C7o4Ac+3XRLrBVNfBJLWD1kj1cdhtwqv/XeNG6gtb3BF62vkf54nqfbJ1Xfuz4AOg9FvqOc1kc7TFxM4erbIx6dzG/bD1ERLA/H97Wk7TEZmbH8l7VFVC485jCsf23vR8FO6Gy6NTfb/WHqCRo3haat4OQFhAUAZHxUF4A342DQa9AfLrj8U76n4dbKMn9bRx63wr44h7vXVfQ+nr7+h7li+vtonX2uj0mvqCi2sZt7y/hl62HCA/y5/1bslRKztTJhluO3ftRUo9x4/BYaHakeDRve/zXEQngV8ev0d4cx+f4dEjIcMbauJeT/bHy1nUFrS949/oe5Yvr7WbrrGLiJiqqbdzxwTIWbD5IaKAf747K5Jw2zc2O1fRKcmHpFOgxqnGN/HTDLYU7wV596ucIjPhtj0eztsd/3awNBGqqfxERV1ExcQOVNTb+/OFy5m88QEiAH1NuzqRHuxZmx3KNo6emdR5QdzE5Otxy3B6Po1/vgMriU79G7XBLuxP3eDRvDyHNm2ZSmIg4x9ist+32PRlfWlfQ+voKX1xvN1hnHWNisqoaO3/+cDmz1+URHGDlnZszueCsGLNjuc7eHHirN9zwKQSGndlwy8n2eDRvB5EJ3newmoiIh9ExJh6g2mbnnn//yux1eQT6W/nHTT5SSo4eaLV7KSx42bHso2tO/T2BEXXs8WjnGG4JCGni0CIi4goeU0y87XThGpud+6bm8O2aXAL9rLx1Y3d6dfSBUgKOMvLL63XfH58O3YYeU0DaNd1wi4iIuBUN5ZjAZjcY80kOM3P2EuBn4c0/dOeSrrFmx2p6NVXwSzbMe/7IrIIWaP872PYDXDkZEs5xPM4bT8cTEfFxGspxUza7wYOfrWBmzl78rRaybzjXN0rJlrnw9YNwcJPjdtJ5MPBFx+m8b/V2lBJvPh1PRETqRcXEhex2g3HTVjJt+R78rBZevf4cLuvm5XsGivbA9+NhzXTH7bCWcOnTkH6dY2jm6FwfIiIiqJi4jGEYPDJzNZ8s3Y3VApNHZDAgNd7sWE2npspxHMkPf4PqMrBYIfN26PswhDT77XFucGqaiIi4DxUTFzAMg8dnreGjRTuxWGDStRkMSk8wO1bT2TrPMWyTv9FxO6knXPEixKed+NiIOJdeg0FERNybikkTMwyDp75cy/sLd2CxwAvD0xlyTmuzYzWN/x22CY2By552XKXSajU3m4iIeAQVkyZkGAYTv1nPlAXbAXhuWCrDuyeaG6op1FTBojccZ9ucathGRETkNDymmHjaPCaGYfDCdxt4a/5WAJ4ZksKIzDYmp2oCW384MmyzwXH7VMM2IiIip6F5TJrIpP9s5JU5jlNjnxzcjZEXtDM3kLMV74XvxsOaaY7boTFw6VOQfr2GbURE5ASax8REr8zZVFtKHr0y2btKia0afnkDfngeqkqPDNvcdmTYxgeuhiwiIk1KxcTJXp+3mUn/cZyNMm5AF27t1d7kRE60bb5j2ObAesftxCwY+HcN24iIiNOomDjR2/O38rdvHcdaPNi/M3/sfZbJiZykeC98/wis/txxW8M2IiLSRFRMnOSdn7Yx4et1ANzfrxOj+55tciInsFXDojdh3nO/Ddv0uBUuHq9hGxERaRIqJk7w/sLtPPXlWgDuvvhs7u3X0eRETnDSYZsXHVf+FRERaSIqJmfoo0U7eWzmGgD+1PssxlzayeREZ6h435Fhm88ct0Ojjwzb3KBhGxERaXIqJmfgkyW7eHj6KgBu69Wev17eGYvFYnKqRrJVw6L/g3kTjxm2uQUufkTDNiIi4jIqJo30+bLd/HXaSgBuvqAd4wd29dxSsu3HI8M2jmNkSMx0TJKWkGFqLBER8T0eU0zcaebXmTl7eOCzFRgG/OG8Njw+KNkzS0nxPvjPo7DqU8ft0Gjo9yRk/F7DNiIiYgrN/NpAX67cyz3//hW7AddnJTFhSCpWq4eVktphm+egqgSw/DZsE9rC7HQiIuKFNPNrE/h29T7u/TgHuwHXdE/0zFKy/Sf46oHfhm1a93CcbZNwjrm5REREUDGpt+/X5HLXR79isxsMO6c1z12d5lmlpCQXvn8UVn3iuB3SAi59EjL+oGEbERFxGyom9fDf9XmM/mg5NXaDwekJvHBNOn6eUkps1bD4LZg7UcM2IiLi9lRMTuOHjQf40wfLqbYZDEyNZ9K1HlRKti+Arx+A/Y7J32jd3XG2Tetzzc0lIiJSBxWTU1iwOZ873l9Klc1O/26xTL4uA38/Dxj2KMmF/zwGK6c6bmvYRkREPISKSR0WbjnIre8tobLGTr+urXj1+nMJcPdSYqs5Mmzz7DHDNqPg4kc1bCMiIh5BxeSIlbsLmfj1esZd0YWKaju3vLuEimo7fTu3JPv35xLo70alpCQXlk5xlI6IOMeyHT87zrbZ75gen4RzYeDfNWwjIiIeRcXkiGnL97Bw60He/GELP2w4wOFqGxd1jOGNP3QnyN/P7HjHK8mFH56DzgMAi2OStGOHbfo9AefcqGEbERHxOD5dTHYXlFNQVo3FAl+s2AvAN6tyMYC0xCieGNyN4AA3KyXHWvUZLH8PKosBC3S/GS55TMM2IiLisXy6mPR6fu4Jy45Og7tydxGX/P0Htj830LWh6lKS6/gAWP+V4/PCVx2fY7pA/2eg46XmZBMREXESj9nXn52dTXJyMpmZmU57zskjMvCv49Rff6uFySMynPZaZ2zpFHirt+Nj/t+Ovy9/Pexeak4uERERJ/L5a+Ws3lPEla/+dMLyL+/uRUrrqDN+fqcpyYXF/wc/TgbsjmX9J0LbCxxfR8T9diCsiIiIm9G1chrIYgHD+O2zW7Hb4Zc3YMFkx+1OA2DjN45SkpBhZjIRERGn8vliEh0eSMvwIOKbBTMiM4mpS3axr7CC6PBAs6M5VFfAzD/D6s8dt/s8DJ0ucxQTERERL+PzxSQ+KoSfxvYl0M+KxWLhhqw2VNns7nGKcPkh+PgG2LkQrP4w+DXIuN4xrNN7rIZuRETE6/h8MQGOKyEWi8U9SsmhbfDhNXBwEwRFwYgPoENvx30RcdB3nLn5REREmoCKiTvavQw+uhbK8yEyEf7wGbTqanYqERGRJqdi4m7WfQmf3wY1hyEuDW74BCLjzU4lIiLiEiom7uSXN+DbcYABHS+D4VMgKNzsVCIiIi6jYuIO7Db4bjwsesNxu/souOJF8NPmERER36J3PrNVlcO022H9l47b/Z6EC+91TKgiIiLiY1RMzFR6AP59HexZCn6BMPRNSLna7FQiIiKmUTExS/4m+HA4FGyHkOZw3b+h7flmpxIRETGViokZdvzsmDjtcAE0bwe//wxiOpqdSkRExHQqJq62+nOY/iewVUHrHnD9xxDe0uxUIiIibkHFxFUMw3ERvtlPOG53uRKGvQ2BoWamEhERcStWswPUV3Z2NsnJyWRmZpodpeFsNfDl/b+VkvP+DNe+r1IiIiLyPyyGYRhmh2iI4uJioqKiKCoqIjIy0uw4p1dZCp+Ngk3fAxa4/Dk4709mpxIREXGp+r5/ayinKRXvc1zzJncl+IfA1f+ArleanUpERMRtqZg0lby1jqsDF++GsJZw/VRI7G52KhEREbemYtIUts6DqTdCZTFEd4Tffwot2pudSkRExO2pmDhbzkcw626w10CbC+C6DyG0hdmpREREPIKKibMYBvzwPMyb6LidcjVc9ToEBJubS0RExIOomDhDTRV8eR/kfOi43et+uPgxsHrM2dgiIiJuQcXkTFUUOY4n2fYDWPxg4N+hxyizU4mIiHgkFZMzUbjLcebNgXUQEAbXvgcdLzU7lYiIiMdSMWmsfSvgw2uhNBfC4+D3n0B8utmpREREPJqKSWNs/B4+vRmqy6BVMtzwCTRLMjuViIiIx1Mxaail78BXD4Bhgw59HNe8CY4yO5WIiIhXUDGpL7sd5jzpuEIwQMbv4crJ4B9oZioRERGvomJSHzWVMONOWP2543afh6H3Q2CxmJtLRETEy6iYnE75Ifj497DzZ7D6w+DXION6s1OJiIh4JRWTUzm0zXE68MFNEBQFIz6ADr3NTiUiIuK1VEzqsnsZfHQtlOdDVJLjQnytupqdSkRExKtpzvSjSnJh7kTH53VfwrsDHaUkPh1um61SIiIi4gLaY3JUSS788BxUlcHC1wADOl4Gw6dAULjZ6URERHyC9pgcZdgdnxe+ChjQ4xa47t8qJSIiIi7k23tMSnIdHxjw5f2/Lc/6o+PMm/J8iIgzLZ6IiIiv8Zg9JtnZ2SQnJ5OZmem8J106Bd7qDW/1gX05vy1f/H+OZUunOO+1RERE5LQshmEYZodoiOLiYqKioigqKiIyMvLMnqx2jwmOi/J9cQ8MeuW3i/FFxGmPiYiIiBPU9/3bt4dyTlY84tMhIcOUOCIiIr7OY4ZyRERExPupmBwVEQe9x2roRkRExES+PZRzrIg46DvO7BQiIiI+TXtMRERExG2omIiIiIjbUDERERERt6FiIiIiIm5DxURERETchoqJiIiIuA0VExEREXEbKiYiIiLiNlRMRERExG2omIiIiIjb8Lgp6Q3DAByXTxYRERHPcPR9++j7eF08rpiUlJQAkJSUZHISERERaaiSkhKioqLqvN9inK66uBm73c7evXu5+OKLWbp06Skfm5mZyZIlS5y6vLi4mKSkJHbt2kVkZGQj16Lx6sroiuep7/ec7nGnur8h2+Z/l2nbaNvUxaxt46ztcqrH6O+ZfmeaQlNsG8MwKCkpISEhAau17iNJPG6PidVqJTExEX9//9NuLD8/v5M+xhnLIyMjTflhqSujK56nvt9zused6v6GbIO6Hqtt0/jHads493mctV1O9Rj9PdPvTFNoqm1zqj0lR3nswa+jR49u9GOctdwMzsrSmOep7/ec7nGnur8h28Cdtgto2zQmj6uYtW2ctV1O9Rj9PdPvTFMwc9t43FCO2YqLi4mKiqKoqMiUFit107ZxX9o27knbxX358rbx2D0mZgkKCuLxxx8nKCjI7CjyP7Rt3Je2jXvSdnFfvrxttMdERERE3Ib2mIiIiIjbUDERERERt6FiIiIiIm5DxURERETchoqJiIiIuA0VEycbOnQozZs3Z/jw4WZHkWPs2rWLPn36kJycTFpaGp9++qnZkQQoLCykR48eZGRkkJKSwttvv212JPkf5eXltG3blgceeMDsKHJEu3btSEtLIyMjg759+5odx+l0urCTzZs3j5KSEt577z0+++wzs+PIEfv27SMvL4+MjAxyc3Pp3r07GzduJCwszOxoPs1ms1FZWUloaChlZWWkpKSwdOlSoqOjzY4mR4wfP57NmzeTlJTEiy++aHYcwVFMVq9eTXh4uNlRmoT2mDhZnz59iIiIMDuG/I/4+HgyMjIAiIuLIyYmhkOHDpkbSvDz8yM0NBSAyspKDMM47SXRxXU2bdrE+vXrGTBggNlRxIeomBxj/vz5DBo0iISEBCwWCzNmzDjhMdnZ2bRr147g4GB69uzJ4sWLXR/UBzlz2yxbtgybzUZSUlITp/Z+ztguhYWFpKenk5iYyIMPPkhMTIyL0ns3Z2ybBx54gIkTJ7oosW9wxnaxWCz07t2bzMxMPvzwQxcldx0Vk2OUlZWRnp5Odnb2Se+fOnUqY8aM4fHHH2f58uWkp6fTv39/9u/f7+KkvsdZ2+bQoUPcdNNNvPXWW66I7fWcsV2aNWvGihUr2LZtGx999BF5eXmuiu/VznTbzJw5k06dOtGpUydXxvZ6zvid+emnn1i2bBmzZs3i2WefZeXKla6K7xqGnBRgTJ8+/bhlWVlZxujRo2tv22w2IyEhwZg4ceJxj5s7d65x9dVXuyKmT2rstqmoqDAuuugi4/3333dVVJ9yJr8zR915553Gp59+2pQxfVJjts3YsWONxMREo23btkZ0dLQRGRlpPPnkk66M7fWc8TvzwAMPGFOmTGnClK6nPSb1VFVVxbJly+jXr1/tMqvVSr9+/Vi4cKGJyaQ+28YwDG6++WYuvvhibrzxRrOi+pT6bJe8vDxKSkoAKCoqYv78+XTu3NmUvL6kPttm4sSJ7Nq1i+3bt/Piiy9y++2389hjj5kV2SfUZ7uUlZXV/s6Ulpby3//+l27dupmSt6n4mx3AU+Tn52Oz2YiNjT1ueWxsLOvXr6+93a9fP1asWEFZWRmJiYl8+umnnH/++a6O61Pqs20WLFjA1KlTSUtLqx3T/eCDD0hNTXV1XJ9Rn+2yY8cO7rjjjtqDXu+++25tExeo798zca36bJe8vDyGDh0KOM5qu/3228nMzHR51qakYuJks2fPNjuCnESvXr2w2+1mx5D/kZWVRU5Ojtkx5DRuvvlmsyPIER06dGDFihVmx2hSGsqpp5iYGPz8/E44MC8vL4+4uDiTUglo27grbRf3pW3jnrRdHFRM6ikwMJDu3bszZ86c2mV2u505c+ZoqMZk2jbuSdvFfWnbuCdtFwcN5RyjtLSUzZs3197etm0bOTk5tGjRgjZt2jBmzBhGjhxJjx49yMrKYvLkyZSVlTFq1CgTU/sGbRv3pO3ivrRt3JO2Sz2YfFaQW5k7d64BnPAxcuTI2se8+uqrRps2bYzAwEAjKyvL+OWXX8wL7EO0bdyTtov70rZxT9oup6dr5YiIiIjb0DEmIiIi4jZUTERERMRtqJiIiIiI21AxEREREbehYiIiIiJuQ8VERERE3IaKiYiIiLgNFRMRERFxGyomIiIi4jZUTERERMRtqJiIiIiI21AxEREREbfx//cmXoS5c/E/AAAAAElFTkSuQmCC",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.xscale(\"log\")\n",
"plt.yscale(\"log\")\n",
"plt.plot(*zip(*res_local), \"-*\")\n",
"plt.plot(*zip(*res_mg), \"-+\")\n",
"plt.legend(['local', 'mg'])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multigrid implementation for higher order spaces \n",
"\n",
"Another important preconditioner that is available in NGSolve is the multigrid preconditioner."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p= 1 : ndof,nsteps= [(8, 2), (21, 4), (65, 7), (225, 8)]\n",
"p= 2 : ndof,nsteps= [(21, 5), (65, 6), (225, 8), (833, 8)]\n",
"p= 3 : ndof,nsteps= [(40, 9), (133, 12), (481, 12), (1825, 13)]\n",
"p= 4 : ndof,nsteps= [(65, 12), (225, 15), (833, 16), (3201, 16)]\n",
"p= 5 : ndof,nsteps= [(96, 14), (341, 19), (1281, 20), (4961, 20)]\n",
"p= 6 : ndof,nsteps= [(133, 16), (481, 23), (1825, 23), (7105, 23)]\n",
"p= 7 : ndof,nsteps= [(176, 18), (645, 25), (2465, 26), (9633, 26)]\n",
"p= 8 : ndof,nsteps= [(225, 19), (833, 27), (3201, 28), (12545, 28)]\n",
"p= 9 : ndof,nsteps= [(280, 20), (1045, 29), (4033, 30), (15841, 30)]\n"
]
}
],
"source": [
"for p in range(1,10):\n",
" r = SolveProblem(h=0.5, p=p, levels=4, condense=False, \n",
" precond=preconditioners.MultiGrid)\n",
" print (\"p=\", p, \": ndof,nsteps=\", r) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We observe that the number of iterations grows mildly with the degree $p$ while remaining bounded with mesh refinement.\n",
"\n",
"Performing static condensation improves the situation:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p= 1 : ndof,nsteps= [(8, 2), (21, 4), (65, 7), (225, 8)]\n",
"p= 2 : ndof,nsteps= [(21, 5), (65, 6), (225, 8), (833, 8)]\n",
"p= 3 : ndof,nsteps= [(40, 5), (133, 6), (481, 7), (1825, 8)]\n",
"p= 4 : ndof,nsteps= [(65, 5), (225, 6), (833, 7), (3201, 8)]\n",
"p= 5 : ndof,nsteps= [(96, 5), (341, 6), (1281, 7), (4961, 8)]\n",
"p= 6 : ndof,nsteps= [(133, 5), (481, 6), (1825, 7), (7105, 8)]\n",
"p= 7 : ndof,nsteps= [(176, 5), (645, 6), (2465, 7), (9633, 8)]\n",
"p= 8 : ndof,nsteps= [(225, 5), (833, 6), (3201, 7), (12545, 8)]\n",
"p= 9 : ndof,nsteps= [(280, 5), (1045, 6), (4033, 7), (15841, 8)]\n"
]
}
],
"source": [
"for p in range(1,10):\n",
" r = SolveProblem(h=0.5, p=p, levels=4, condense=True, \n",
" precond=preconditioners.MultiGrid)\n",
" print (\"p=\", p, \": ndof,nsteps=\", r) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Element-wise BDDC preconditioner\n",
"\n",
"A built-in element-wise BDDC (Balancing Domain Decomposition preconditioner with Constraints) preconditioner is also available. In contrast to local or multigrid preconditioners, the BDDC preconditioner needs access to the element matrices. This is exactly why we need to register the preconditioner with the bilinear form `bfa` before calling `bfa.Assemble()`."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p= 1 : ndof,nsteps= [(8, 2), (21, 2), (65, 2), (225, 2)]\n",
"p= 2 : ndof,nsteps= [(21, 7), (65, 12), (225, 16), (833, 17)]\n",
"p= 3 : ndof,nsteps= [(40, 11), (133, 16), (481, 20), (1825, 22)]\n",
"p= 4 : ndof,nsteps= [(65, 13), (225, 18), (833, 24), (3201, 26)]\n",
"p= 5 : ndof,nsteps= [(96, 14), (341, 19), (1281, 26), (4961, 29)]\n",
"p= 6 : ndof,nsteps= [(133, 14), (481, 20), (1825, 28), (7105, 31)]\n",
"p= 7 : ndof,nsteps= [(176, 15), (645, 21), (2465, 30), (9633, 33)]\n",
"p= 8 : ndof,nsteps= [(225, 16), (833, 23), (3201, 31), (12545, 35)]\n",
"p= 9 : ndof,nsteps= [(280, 16), (1045, 24), (4033, 33), (15841, 36)]\n"
]
}
],
"source": [
"for p in range(1,10):\n",
" r = SolveProblem(h=0.5, p=p, levels=4, condense=True, \n",
" precond=preconditioners.BDDC)\n",
" print (\"p=\", p, \": ndof,nsteps=\", r) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The BDDC preconditioner needs more iterations, but the work per iteration is less, so performance is similar to multigrid. **This element-wise BDDC preconditioner works well in shared memory parallel as well as in distributed memory mode.** See $\\S$[2.1.4](../unit-2.1.4-bddc/bddc.ipynb) for more about BDDC preconditioner and how to combine it with an algebraic multigrid coarse solver. "
]
}
],
"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.12.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
4426a015fThis file is pinned to the selected release. Download the original to open it in your design tool.