adapMP

notebooks/IGA_FEM.ipynb

@@ -0,0 +1,340 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "8ba4e535-7510-4e12-bf36-f9d5e9ee2ed9",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Populating the interactive namespace from numpy and matplotlib\n"
+     ]
+    }
+   ],
+   "source": [
+    "%pylab inline\n",
+    "import numpy as np\n",
+    "import scipy\n",
+    "import math\n",
+    "import time\n",
+    "from scipy.sparse.linalg import onenormest, splu, LinearOperator\n",
+    "from scipy.sparse import coo_matrix, block_diag, identity, hstack\n",
+    "import matplotlib.pyplot as plt\n",
+    "from pyiga import assemble, bspline, vform, geometry, vis, solvers, utils, topology\n",
+    "#from patchmesh import *\n",
+    "#from sksparse.cholmod import cholesky\n",
+    "#from patchmesh3D import *\n",
+    "#from multipatch import *\n",
+    "\n",
+    "np.set_printoptions(linewidth=100000)\n",
+    "np.set_printoptions(precision=2)\n",
+    "np.set_printoptions(formatter={'float': lambda x: \"{0: 0.3f}\".format(x) if x>=0 else \"{0:0.3f}\".format(x)})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "4c225d79-ee4e-4581-80b1-0927872fc136",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def condest(A):\n",
+    "    luA = splu(A)\n",
+    "    iA = LinearOperator(luA.shape, matvec = lambda x : luA.solve(x), rmatvec = lambda x : luA.solve(x))\n",
+    "    return onenormest(iA)*onenormest(A)\n",
+    "\n",
+    "def pcg(Afuns, f, tol = 1e-5, maxit = 100, pfuns = 1, output = False):     \n",
+    "    maxit = int(maxit)\n",
+    "    if not callable(pfuns):\n",
+    "        pfun = lambda x : x\n",
+    "        # splu_pfun = sp.linalg.splu(pfuns,permc_spec='COLAMD')\n",
+    "        # pfun = lambda x : splu_pfun.solve(x)\n",
+    "    else:\n",
+    "        pfun = pfuns\n",
+    "    if not callable(Afuns):\n",
+    "        Afun = lambda x : Afuns@x\n",
+    "    else:\n",
+    "        Afun = Afuns\n",
+    "    if not isinstance(f,np.ndarray):\n",
+    "        d = f.A.ravel()\n",
+    "    else:\n",
+    "        d = f.ravel() \n",
+    "    # print('Cond about',condest(pfuns@Afuns))\n",
+    "    w = pfun(d)\n",
+    "    rho = w@d\n",
+    "    err0 = np.sqrt(rho)\n",
+    "    s = w\n",
+    "    u = 0*d.copy()\n",
+    "    for it in range(maxit):\n",
+    "        As = Afun(s)\n",
+    "        alpha = rho/(As@s)\n",
+    "        u = u + alpha*s\n",
+    "        d = d - alpha*As\n",
+    "        w = pfun(d)\n",
+    "        rho1 = rho\n",
+    "        rho = w@d\n",
+    "        err = np.sqrt(rho)\n",
+    "        if err < tol*err0:\n",
+    "            break\n",
+    "        beta = rho/rho1\n",
+    "        s = w + beta*s\n",
+    "    if output:\n",
+    "        print('pcg stopped after ' + str(it) + ' iterations with relres ' + str(err/err0))\n",
+    "    return u,it,d"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "bc18c467-f277-4911-83bb-b72191d569e0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def refine(kvs, mult=1):\n",
+    "    return tuple(kv.h_refine(mult=mult) for kv in kvs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 33,
+   "id": "5818aac7-8be0-4c6c-a359-01f4c9085b08",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "u = lambda x,y: exp(x**2+y**2)\n",
+    "f = lambda x,y: -4*(x**2+y**2+1)*exp(x**2+y**2)\n",
+    "geo=geometry.unit_square()\n",
+    "maxiter=8\n",
+    "degree = (1,2,3,4)\n",
+    "numdofs, assembly_t, nnz, nnzq, nnzL, nnzLq, cond, cg_it = np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter)\n",
+    "\n",
+    "for p in degree:\n",
+    "    kvs = 2*(bspline.make_knots(p,0.0,1.0,2),)\n",
+    "    for i in range(maxiter):\n",
+    "        numdofs[(p-1)*maxiter+i]=np.prod([kv.numdofs for kv in kvs])\n",
+    "        t=time.time()\n",
+    "        Ah = assemble.stiffness(kvs, geo)\n",
+    "        Fh = assemble.inner_products(kvs, f, f_physical=True, geo=geo).ravel()\n",
+    "        assembly_t[(p-1)*maxiter+i]=time.time()-t\n",
+    "        nnz[(p-1)*maxiter+i]=Ah.nnz#/np.prod(Ah.shape)\n",
+    "        nnzq[(p-1)*maxiter+i]=100*Ah.nnz/np.prod(Ah.shape)\n",
+    "        \n",
+    "        bcs = assemble.compute_dirichlet_bcs(kvs, geo, [('bottom', u), ('top', u), ('left', u), ('right', u)])\n",
+    "        LS = assemble.RestrictedLinearSystem(Ah, Fh, bcs)\n",
+    "        solver = scipy.sparse.linalg.splu(LS.A)\n",
+    "        nnzL[(p-1)*maxiter+i] = solver.L.nnz#/np.prod(solver.L.shape)\n",
+    "        nnzLq[(p-1)*maxiter+i] = 100*solver.L.nnz/np.prod(solver.L.shape)\n",
+    "        cond[(p-1)*maxiter+i] = condest(LS.A)\n",
+    "        _, it, _ = pcg(LS.A,LS.b,maxit=1000)\n",
+    "        cg_it[(p-1)*maxiter+i] = it\n",
+    "        #bcs = assemble.compute_dirichlet_bcs(kvs, geo, [('bottom',0),('right',0),('left', 0), ('top', 0)])\n",
+    "        if i!=maxiter-1:\n",
+    "            kvs = refine(kvs)\n",
+    "np.savetxt('Iga.txt',np.c_[numdofs, assembly_t, nnz, nnzq, nnzL, nnzLq, cond, cg_it], fmt=('%d', '%1.3f','%d','%1.3f','%d','%1.3f', '%1.3e','%d'), delimiter = \" & \")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "id": "8366adbe-4051-4c5a-8ff4-2ac92a71fbcc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "numdofs, assembly_t, nnz, nnzL, cond, cg_it = np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter), np.zeros(len(degree)*maxiter)\n",
+    "\n",
+    "for p in degree:\n",
+    "    kvs = 2*(bspline.make_knots(p,0.0,1.0,2,mult=p),)\n",
+    "    for i in range(maxiter):\n",
+    "        numdofs[(p-1)*maxiter+i]=np.prod([kv.numdofs for kv in kvs])\n",
+    "        t=time.time()\n",
+    "        Ah = assemble.stiffness(kvs, geo)\n",
+    "        Fh = assemble.inner_products(kvs, f, f_physical=True, geo=geo).ravel()\n",
+    "        assembly_t[(p-1)*maxiter+i]=time.time()-t\n",
+    "        nnz[(p-1)*maxiter+i]=Ah.nnz#/np.prod(Ah.shape)\n",
+    "        nnzq[(p-1)*maxiter+i]=100*Ah.nnz/np.prod(Ah.shape)\n",
+    "        \n",
+    "        bcs = assemble.compute_dirichlet_bcs(kvs, geo, [('bottom', u), ('top', u), ('left', u), ('right', u)])\n",
+    "        LS = assemble.RestrictedLinearSystem(Ah, Fh, bcs)\n",
+    "        solver = scipy.sparse.linalg.splu(LS.A)\n",
+    "        nnzL[(p-1)*maxiter+i] = solver.L.nnz#/np.prod(solver.L.shape)\n",
+    "        nnzLq[(p-1)*maxiter+i] = 100*solver.L.nnz/np.prod(solver.L.shape)\n",
+    "        cond[(p-1)*maxiter+i] = condest(LS.A)\n",
+    "        _, it, _ = pcg(LS.A,LS.b,maxit=1000)\n",
+    "        cg_it[(p-1)*maxiter+i] = it\n",
+    "        #bcs = assemble.compute_dirichlet_bcs(kvs, geo, [('bottom',0),('right',0),('left', 0), ('top', 0)])\n",
+    "        if i!=maxiter-1:\n",
+    "            kvs = refine(kvs,mult=p)\n",
+    "np.savetxt('FEM.txt',np.c_[numdofs, assembly_t, nnz, nnzq, nnzL, nnzLq, cond, cg_it], fmt=('%d', '%1.3f','%d','%1.3f','%d','%1.3f', '%1.3e','%d'), delimiter = \" & \")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "8202037e-21ed-454b-af0d-bd3865d2c8eb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Ah = assemble.stiffness(2*(bspline.make_knots(4,0.0,1.0,10,mult=1),), geo=geometry.unit_cube(2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "id": "0a2037db-d24b-43ba-a344-d5bb80bfc058",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "3e+00\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(\"%1.0e\" % 300.0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "cdad44ce-b26d-4212-92a5-061907ffdd81",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1191016\n",
+      "0.15817920254315804\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(Ah.nnz)\n",
+    "print(Ah.nnz/(np.prod(Ah.shape)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "7fcfe66e-b2b6-411d-b316-87d16d7033a4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Ah = assemble.stiffness(3*(bspline.make_knots(4,0.0,1.0,10,mult=4),), geo=geometry.unit_cube(3))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "bda47b3d-6cc3-4a32-a7b1-431e45a74fe1",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "13997521\n",
+      "0.0029467818577920764\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(Ah.nnz)\n",
+    "print(Ah.nnz/(np.prod(Ah.shape)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 181,
+   "id": "22410362-ccc5-42cf-9a13-0678d36bf407",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kvs = 2*(bspline.make_knots(2,0.0,1.0,3,mult=2),)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 182,
+   "id": "ed4e1009-6731-406e-97e8-82eb5e0674ab",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "4"
+      ]
+     },
+     "execution_count": 182,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "kvs[0].first_active_at(0.75)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 184,
+   "id": "837f8a5d-12c2-4ca7-9e8f-a51f51f88e1a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "7"
+      ]
+     },
+     "execution_count": 184,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "kvs[0].numdofs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "64e003d0-e5a1-4083-8fc8-c85b6564219c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.sparse.linalg import onenormest, splu, LinearOperator\n",
+    "\n",
+    "def condest(A):\n",
+    "    luA = splu(A)\n",
+    "    iA = LinearOperator(luA.shape, matvec = lambda x : luA.solve(x), rmatvec = lambda x : luA.solve(x))\n",
+    "    return onenormest(iA)*onenormest(A)"
+   ]
+  }
+ ],
+ "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.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}