demo 3d.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting demo_3d.py\n"
]
}
],
"source": [
"%%file demo_3d.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",
" rhs[..., 1] = -9.81# + rand\n",
"\n",
"OptDB = PETSc.Options()\n",
"Lx = OptDB.getInt('Lx', 10)\n",
"Ly = OptDB.getInt('Ly', 1)\n",
"Lz = OptDB.getInt('Lz', 1)\n",
"n = OptDB.getInt('n', 16)\n",
"nx = OptDB.getInt('nx', Lx*n)\n",
"ny = OptDB.getInt('ny', Ly*n)\n",
"nz = OptDB.getInt('nz', Lz*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",
"hz = Lz/(nz - 1)\n",
"h = [hx, hy, hz]\n",
"\n",
"da = PETSc.DMDA().create([nx, ny, nz], dof=3, stencil_width=1)\n",
"da.setUniformCoordinates(xmax=Lx, ymax=Ly, zmax=Lz)\n",
"da.setMatType(PETSc.Mat.Type.IS)\n",
"\n",
"path = './output_3d/'\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",
"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",
"def lame_coeff(x, y, z, v1, v2):\n",
" output = np.empty(x.shape)\n",
" mask = np.logical_or(np.logical_and(.2<=z, z<=.4),np.logical_and(.6<=z, z<=.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",
"x = da.createGlobalVec()\n",
"b = buildRHS(da, h, rhs)\n",
"A = buildElasticityMatrix(da, h, lamb, mu)\n",
"A.assemble()\n",
"\n",
"bcApplyWest(da, A, b)\n",
"bcopy = b.copy()\n",
"\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.setInitialGuessNonzero(True)\n",
"ksp.setFromOptions()\n",
"\n",
"ksp.solve(b, x)\n",
"\n",
"viewer = PETSc.Viewer().createVTK(path + 'solution_3d.vts', 'w', comm = PETSc.COMM_WORLD)\n",
"x.view(viewer)\n",
"\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], coords_a[:, :, :, 2], E1, E2)\n",
"nu = lame_coeff(coords_a[:, :, :, 0], coords_a[:, :, :, 1], coords_a[:, :, :, 2], 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",
"viewer.destroy()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Subdomain number 3 contributes 6 coarse vectors as zero energy modes of local solver\n",
"Subdomain number 3 contributes 6 coarse vectors in total\n",
"Subdomain number 2 contributes 6 coarse vectors as zero energy modes of local solver\n",
"Subdomain number 2 contributes 6 coarse vectors in total\n",
"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 6 coarse vectors as zero energy modes of local solver\n",
"Subdomain number 1 contributes 6 coarse vectors in total\n",
"There are 18 vectors in the coarse space.\n",
"multipreconditioning initial iteration\n",
" 0 KSP Residual norm 1.607932673323e+02 \n",
"\tnatural_norm -> 5.59670982e+00\n",
"\tti -> 0.00000000e+00\n",
"multipreconditioning this iteration\n",
" 1 KSP Residual norm 1.292233160585e+02 \n",
"\tnatural_norm -> 3.66856270e+01\n",
"\tti -> 1.08805060e-02\n",
"multipreconditioning this iteration\n",
" 2 KSP Residual norm 1.412616300838e+02 \n",
"\tnatural_norm -> 7.91353302e+01\n",
"\tti -> 1.85911544e-05\n",
"multipreconditioning this iteration\n",
" 3 KSP Residual norm 4.380698167085e+02 \n",
"\tnatural_norm -> 2.48213524e+02\n",
"\tti -> 6.74592702e-05\n",
"multipreconditioning this iteration\n",
" 4 KSP Residual norm 4.446605915798e+02 \n",
"\tnatural_norm -> 9.22485277e+01\n",
"\tti -> 1.76134061e-04\n",
"multipreconditioning this iteration\n",
" 5 KSP Residual norm 1.530281901669e+02 \n",
"\tnatural_norm -> 2.64051339e+01\n",
"\tti -> 2.27470131e-03\n",
"multipreconditioning this iteration\n",
" 6 KSP Residual norm 3.053836608546e+02 \n",
"\tnatural_norm -> 4.24387313e+01\n",
"\tti -> 1.35905736e-03\n",
" 7 KSP Residual norm 2.025618136576e+01 \n",
"\tnatural_norm -> 2.86136488e-01\n",
"\tti -> 7.45527474e+01\n",
" 8 KSP Residual norm 1.633561158521e+01 \n",
"\tnatural_norm -> 1.62756942e-01\n",
"\tti -> 5.68544944e-01\n",
" 9 KSP Residual norm 5.192570025925e+00 \n",
"\tnatural_norm -> 8.93538483e-02\n",
"\tti -> 5.13398502e-01\n",
" 10 KSP Residual norm 3.804896305678e+00 \n",
"\tnatural_norm -> 5.70709200e-02\n",
"\tti -> 5.16455326e-01\n",
" 11 KSP Residual norm 2.569423123279e+00 \n",
"\tnatural_norm -> 3.14418164e-02\n",
"\tti -> 4.46254988e-01\n",
" 12 KSP Residual norm 2.334275734190e+00 \n",
"\tnatural_norm -> 2.37510957e-02\n",
"\tti -> 5.91126718e-01\n",
" 13 KSP Residual norm 5.529487252972e-01 \n",
"\tnatural_norm -> 9.28772018e-03\n",
"\tti -> 1.28485150e+00\n",
" 14 KSP Residual norm 3.986599842461e-01 \n",
"\tnatural_norm -> 6.26572422e-03\n",
"\tti -> 3.72583445e-01\n",
" 15 KSP Residual norm 1.747913265366e-01 \n",
"\tnatural_norm -> 1.73984334e-03\n",
"\tti -> 1.75735096e+00\n",
" 16 KSP Residual norm 1.101212440906e-01 \n",
"\tnatural_norm -> 1.28333379e-03\n",
"\tti -> 5.64881232e-01\n",
" 17 KSP Residual norm 7.729160714217e-02 \n",
"\tnatural_norm -> 8.52182649e-04\n",
"\tti -> 3.42558277e-01\n",
" 18 KSP Residual norm 3.139783442015e-02 \n",
"\tnatural_norm -> 2.89672134e-04\n",
"\tti -> 1.47552733e+00\n",
" 19 KSP Residual norm 1.880284157141e-02 \n",
"\tnatural_norm -> 1.65763072e-04\n",
"\tti -> 6.18623758e-01\n",
" 20 KSP Residual norm 1.311832777508e-02 \n",
"\tnatural_norm -> 1.14162940e-04\n",
"\tti -> 4.93348725e-01\n",
" 21 KSP Residual norm 5.486272159507e-03 \n",
"\tnatural_norm -> 4.34202281e-05\n",
"\tti -> 9.31707361e-01\n",
" 22 KSP Residual norm 3.348948267277e-03 \n",
"\tnatural_norm -> 3.10264661e-05\n",
"\tti -> 4.71945626e-01\n",
" 23 KSP Residual norm 1.997199094256e-03 \n",
"\tnatural_norm -> 1.91982388e-05\n",
"\tti -> 4.67892538e-01\n",
" 24 KSP Residual norm 6.759517088286e-04 \n",
"\tnatural_norm -> 6.65257265e-06\n",
"\tti -> 1.43222695e+00\n",
" 25 KSP Residual norm 5.090408497775e-04 \n",
"\tnatural_norm -> 4.98139457e-06\n",
"\tti -> 2.58189484e-01\n",
" 26 KSP Residual norm 3.388640771041e-04 \n",
"\tnatural_norm -> 2.37443890e-06\n",
"\tti -> 1.20625254e+00\n",
" 27 KSP Residual norm 1.222987667731e-04 \n",
"\tnatural_norm -> 8.67066551e-07\n",
"\tti -> 1.81789157e+00\n",
" 28 KSP Residual norm 6.531347718113e-05 \n",
"\tnatural_norm -> 6.73054653e-07\n",
"\tti -> 4.73069334e-01\n",
" 29 KSP Residual norm 5.287930795043e-05 \n",
"\tnatural_norm -> 4.44224709e-07\n",
"\tti -> 3.30004865e-01\n",
" 30 KSP Residual norm 1.923467995097e-05 \n",
"\tnatural_norm -> 1.55135222e-07\n",
"\tti -> 1.62217244e+00\n",
" 31 KSP Residual norm 1.277966311929e-05 \n",
"\tnatural_norm -> 8.84535505e-08\n",
"\tti -> 9.26566133e-01\n",
" 32 KSP Residual norm 6.672587978298e-06 \n",
"\tnatural_norm -> 4.34644849e-08\n",
"\tti -> 9.05676258e-01\n",
" 33 KSP Residual norm 3.622973713319e-06 \n",
"\tnatural_norm -> 2.70984183e-08\n",
"\tti -> 5.68387245e-01\n"
]
}
],
"source": [
"!mpiexec -np 4 python demo_3d.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": "ef6f175ec2c84ca3bcd0278c43fbbfad",
"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_3d', 'solution_3d.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
}
