{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Overwriting demo_2d.py\n" ] } ], "source": [ "%%file demo_2d.py\n", "\n", "from __future__ import print_function, division\n", "import os\n", "import sys, petsc4py\n", "petsc4py.init(sys.argv)\n", "import mpi4py.MPI as mpi\n", "from petsc4py import PETSc\n", "import numpy as np\n", "from GenEO import *\n", "\n", "def rhs(coords, rhs):\n", " n = rhs.shape\n", " rhs[..., 1] = -9.81\n", "\n", "OptDB = PETSc.Options()\n", "Lx = OptDB.getInt('Lx', 10)\n", "Ly = OptDB.getInt('Ly', 1)\n", "n = OptDB.getInt('n', 16)\n", "nx = OptDB.getInt('nx', Lx*n)\n", "ny = OptDB.getInt('ny', Ly*n)\n", "E1 = OptDB.getReal('E1', 10**6)\n", "E2 = OptDB.getReal('E2', 1)\n", "nu1 = OptDB.getReal('nu1', 0.4)\n", "nu2 = OptDB.getReal('nu2', 0.4)\n", "\n", "hx = Lx/(nx - 1)\n", "hy = Ly/(ny - 1)\n", "\n", "da = PETSc.DMDA().create([nx, ny], dof=2, stencil_width=1)\n", "da.setUniformCoordinates(xmax=Lx, ymax=Ly)\n", "da.setMatType(PETSc.Mat.Type.IS)\n", "da.setFieldName(0, 'u')\n", "da.setFieldName(1, 'v')\n", "\n", "path = './output_2d/'\n", "if mpi.COMM_WORLD.rank == 0:\n", " if not os.path.exists(path):\n", " os.mkdir(path)\n", " else:\n", " os.system('rm {}/*.vts'.format(path))\n", "\n", "def lame_coeff(x, y, v1, v2):\n", " output = np.empty(x.shape)\n", " mask = np.logical_or(np.logical_and(.2<=y, y<=.4),np.logical_and(.6<=y, y<=.8))\n", " output[mask] = v1\n", " output[np.logical_not(mask)] = v2\n", " return output\n", "\n", "# non constant Young's modulus and Poisson's ratio \n", "E = buildCellArrayWithFunction(da, lame_coeff, (E1,E2))\n", "nu = buildCellArrayWithFunction(da, lame_coeff, (nu1,nu2))\n", "\n", "lamb = (nu*E)/((1+nu)*(1-2*nu)) \n", "mu = .5*E/(1+nu)\n", "\n", "class callback:\n", " def __init__(self, da):\n", " self.da = da\n", " ranges = da.getRanges()\n", " ghost_ranges = da.getGhostRanges()\n", " \n", " self.slices = []\n", " for r, gr in zip(ranges, ghost_ranges):\n", " self.slices.append(slice(gr[0], r[1]))\n", " self.slices = tuple(self.slices)\n", "\n", " self.it = 0\n", "\n", " def __call__(self, locals):\n", " pyKSP = locals['self']\n", " proj = pyKSP.mpc.proj\n", "\n", " viewer_x = PETSc.Viewer().createVTK(path + 'cg_ite_{}.vts'.format(self.it), 'w', comm = PETSc.COMM_WORLD)\n", " locals['x'].view(viewer_x)\n", " viewer_x.destroy()\n", " \n", " if self.it == 0:\n", " work, _ = proj.A.getVecs()\n", " for i, vec in enumerate(proj.coarse_vecs):\n", " if vec:\n", " proj.workl = vec.copy()\n", " else:\n", " proj.workl.set(0.)\n", " work.set(0)\n", " proj.scatter_l2g(proj.workl, work, PETSc.InsertMode.ADD_VALUES)\n", "\n", " viewer = PETSc.Viewer().createVTK(path + 'coarse_vec_{}.vts'.format(i), 'w', comm = PETSc.COMM_WORLD)\n", " tmp = self.da.createGlobalVec()\n", " tmpl_a = self.da.getVecArray(tmp)\n", " work_a = self.da.getVecArray(work)\n", " tmpl_a[:] = work_a[:]\n", " tmp.view(viewer)\n", " viewer.destroy()\n", " self.it += 1\n", "\n", "\n", "x = da.createGlobalVec()\n", "b = buildRHS(da, [hx, hy], rhs)\n", "A = buildElasticityMatrix(da, [hx, hy], lamb, mu)\n", "A.assemble()\n", "bcApplyWest(da, A, b)\n", "\n", "#Setup the preconditioner (or multipreconditioner) and the coarse space\n", "pcbnn = PCBNN(A)\n", "\n", "# Set initial guess\n", "x.setRandom()\n", "xnorm = b.dot(x)/x.dot(A*x)\n", "x *= xnorm\n", "\n", "ksp = PETSc.KSP().create()\n", "ksp.setOperators(A)\n", "ksp.setType(ksp.Type.PYTHON)\n", "pyKSP = KSP_AMPCG(pcbnn)\n", "pyKSP.callback = callback(da)\n", "ksp.setPythonContext(pyKSP)\n", "ksp.setFromOptions()\n", "ksp.setInitialGuessNonzero(True)\n", "\n", "ksp.solve(b, x)\n", "\n", "def write_simu_info(da, viewer):\n", " lamb_petsc = da.createGlobalVec()\n", " lamb_a = da.getVecArray(lamb_petsc)\n", " coords = da.getCoordinates()\n", " coords_a = da.getVecArray(coords)\n", " E = lame_coeff(coords_a[:, :, 0], coords_a[:, :, 1], E1, E2)\n", " nu = lame_coeff(coords_a[:, :, 0], coords_a[:, :, 1], nu1, nu2)\n", "\n", " lamb_a[:, :, 0] = (nu*E)/((1+nu)*(1-2*nu)) \n", " lamb_a[:, :, 1] = mpi.COMM_WORLD.rank\n", " lamb_petsc.view(viewer)\n", "\n", "\n", "viewer = PETSc.Viewer().createVTK(path + 'solution_2d.vts', 'w', comm = PETSc.COMM_WORLD)\n", "x.view(viewer)\n", "write_simu_info(da, viewer)\n", "viewer.destroy()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Subdomain number 0 contributes 0 coarse vectors as zero energy modes of local solver\n", "Subdomain number 0 contributes 0 coarse vectors in total\n", "Subdomain number 1 contributes 3 coarse vectors as zero energy modes of local solver\n", "Subdomain number 1 contributes 3 coarse vectors in total\n", "Subdomain number 3 contributes 3 coarse vectors as zero energy modes of local solver\n", "Subdomain number 3 contributes 3 coarse vectors in total\n", "Subdomain number 2 contributes 3 coarse vectors as zero energy modes of local solver\n", "Subdomain number 2 contributes 3 coarse vectors in total\n", "There are 9 vectors in the coarse space.\n", "multipreconditioning initial iteration\n", " 0 KSP Residual norm 8.630197784840e+02 \n", "\tnatural_norm -> 2.00462014e+02\n", "\tti -> 0.00000000e+00\n", "multipreconditioning this iteration\n", " 1 KSP Residual norm 2.865376542965e+03 \n", "\tnatural_norm -> 4.37483489e+01\n", "\tti -> 7.47459199e-02\n", "multipreconditioning this iteration\n", " 2 KSP Residual norm 2.900386632548e+03 \n", "\tnatural_norm -> 2.95866047e+01\n", "\tti -> 2.30814813e-04\n", " 3 KSP Residual norm 1.977093599955e+03 \n", "\tnatural_norm -> 1.02475612e+01\n", "\tti -> 5.65467225e-01\n", " 4 KSP Residual norm 8.838191282536e+02 \n", "\tnatural_norm -> 9.37288194e+00\n", "\tti -> 4.89892005e-01\n", " 5 KSP Residual norm 3.865905875349e+02 \n", "\tnatural_norm -> 2.35257625e+00\n", "\tti -> 3.99665534e+00\n", " 6 KSP Residual norm 8.724468700199e+01 \n", "\tnatural_norm -> 7.43514246e-01\n", "\tti -> 5.03562243e+00\n", " 7 KSP Residual norm 6.672809835644e+01 \n", "\tnatural_norm -> 2.83624035e-01\n", "\tti -> 1.11134346e+00\n", " 8 KSP Residual norm 1.886871132444e+01 \n", "\tnatural_norm -> 5.34966961e-02\n", "\tti -> 1.18963143e+01\n", " 9 KSP Residual norm 3.518634563815e+00 \n", "\tnatural_norm -> 2.65417910e-02\n", "\tti -> 2.86041402e+00\n", " 10 KSP Residual norm 2.431257845448e+00 \n", "\tnatural_norm -> 1.55113887e-02\n", "\tti -> 6.35645347e-01\n", " 11 KSP Residual norm 1.336323965736e+00 \n", "\tnatural_norm -> 1.67404683e-03\n", "\tti -> 2.60090408e+01\n", " 12 KSP Residual norm 3.661523265534e-01 \n", "\tnatural_norm -> 3.86927633e-04\n", "\tti -> 9.40814819e+00\n", " 13 KSP Residual norm 2.092563511534e-02 \n", "\tnatural_norm -> 7.11242725e-05\n", "\tti -> 1.70597722e+01\n", " 14 KSP Residual norm 1.047012961882e-02 \n", "\tnatural_norm -> 1.70066913e-05\n", "\tti -> 8.37961941e+00\n", " 15 KSP Residual norm 5.351566510790e-03 \n", "\tnatural_norm -> 1.56932469e-05\n", "\tti -> 4.26618612e-01\n", " 16 KSP Residual norm 1.084615874674e-03 \n", "\tnatural_norm -> 3.39244216e-06\n", "\tti -> 5.53096481e+00\n", " 17 KSP Residual norm 4.625807632615e-04 \n", "\tnatural_norm -> 5.65141702e-07\n", "\tti -> 2.05101245e+01\n", " 18 KSP Residual norm 7.948087158020e-05 \n", "\tnatural_norm -> 2.78209343e-07\n", "\tti -> 2.41692183e+00\n", "multipreconditioning this iteration\n", " 19 KSP Residual norm 6.744891661848e-06 \n", "\tnatural_norm -> 2.69757887e-07\n", "\tti -> 4.97777267e-02\n" ] } ], "source": [ "!mpiexec -np 4 python demo_2d.py -AMPCG_verbose -ksp_monitor -PCBNN_verbose" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from plot import plot_solution" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "ba261530d25b4a4b915d5560b5c26f43", "version_major": 2, "version_minor": 0 }, "text/plain": [ "HBox(children=(Renderer(camera=PerspectiveCamera(position=(0.0, 0.0, 13.0), quaternion=(0.0, 0.0, 0.0, 1.0), s…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_solution('output_2d', 'solution_2d.vts')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }