demo 2d.ipynb
{
"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
}
